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ABSTRACT 


?  Two  papers  given  4>y—$ACLANTCEN  personnel}  at  the  NATO  Advanced  Research 
Workshop  on  Hybrid  Formulation  of  Wave  Propagation  and  Scatteri ng  4n  Castel 
^—c  fiandol fcr,  Romei  fttrf-y  on  30  Aug€fet}to  3  September^  1983  are  published  as  a 
Review  of  SACLANTCEN's  recent  programme  in  underwater-acoustic  modelling. 
The  first  briefly  reviews  the  physics  of  sound  propagation  in  the  ocean. 
In  it  the  mathematical  foundations  of  the  most  widely  used  acoustic  models 
(ray,  mode,  fast  field,  parabolic  equation)  are  presented  and  the  areas  of 
applicability  of  the  various  models  are  indicated.  A  few  numerical 

examples  are  Included  to  show  the  consistency  among  the  different  computer 
models  in  overlapping  regimes  of  validity.  A  series  of  computational 
examples  is  given  to  demonstrate  the  applicability  of  these  models  to  a 
wide  range  of  general  wave-propagation  problems.  The  second  paper  presents 
a  new  numerical  model,  of  the  fast  field  type,  where  the  depth-separated 
wave  equation  is  solved  by  a  numerical  technique  very  similar  to  that  used 
in  finite-element  programs.  The  speed  Improvement  over  existing  models  of 
the  same  type  is  considerable,  especially  In  cases  with  many  sources  and 
receivers.  The  model  has  been  used  for  studying  both  seismic  pulse  propa¬ 
gation  in  shallow  water  and  reflection  of  pulsed  ultrasonic  beams  from  a 
fluid/solid  interface. 
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PREFACE 


This  report  contains  two  SACLANTCEN  contributions  to  a  NATO  Advanced 
Research  Workshop  on  Hybrid  Formulation  of  Wave  Propagation  and  Scattering, 
held  in  Castel  Gandolfo,  Rome,  Italy,  on  30  August  to  3  September  1983. 
The  Workshop  was  directed  by  Professor  L.B.  Felsen,  Polytechnic  Institute 
of  New  York,  with  the  scope  of  bringing  together  researchers  in  different 
fields  of  wave  propagation  (electromagnetics,  optics,  seismics,  underwater 
acoustics)  to  exchange  ideas  and  approaches  to  solving  complex  wave- 
propagation  problems. 

A. 

SACLANTCEN  was  invited  to  present  two  papers:  a  review  of  numerical  models 
in  underwater  acoustics  (F.B.  Jensen),  and  a  recently  developed  solution 
technique  of  the  fast-field  type  (H.  Schmidt).  These  two  papers  are  here 
published  as  a  single  report,  since  it  is  felt  that  together  they  give  a 
good  representation  of  tne  successful  research  programme  of  SACLANTCEN 's 
Environmental  Modelling  Group  over  the  past  years. 

The  full  proceedings  of  the  Workshop,  containing  24  contributions,  is 
available  in  <1>. 


Scattering. 


ed.  Hybrid  Formulation  of  Wave  Propagation 
Dordrecht,  Netherlands,  Martinus  Nijhoff,  1984. 


and 
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NUMERICAL  MODELS  IN  UNDERWATER  ACOUSTICS 


Finn  B.  Jensen 

SACLANT  ASW  Research  Centre 
Vlale  San  Bartolomeo  400 
19026  La  Spezia,  Italy 


ABSTRACT 

The  physics  of  sound  propagation  In  the  ocean  Is  briefly 
reviewed.  The  mathematical  foundation  of  the  most  widely  used 
acoustic  models  (ray,  mode,  fast  field,  parabolic  equation)  Is 
presented  and  the  areas  of  applicability  of  the  various  models  are 
indicated.  A  few  numerical  examples  are  Included  to  show  the  con¬ 
sistency  among  the  different  computer  models  in  overlapping  regi¬ 
mes  of  validity.  Finally,  we  show  a  series  of  computational 
examples  that  demonstrate  the  applicability  of  these  models  to  a 
wide  range  of  general  wave-propagation  problems. 


INTRODUCTION 

The  modern  era  of  underwater  acoustics  essentially  dates  back 
to  the  beginning  of  World  War  II,  where  considerable  effort  went 
into  Improving  submarine  detection  by  acoustic  means.  This  effort 
has  continued,  promoted  by  naval  Interests  in  developing  still 
better  and  more  reliable  sonar  systems.  To  achieve  optimum  sonar 
design  one  needs  to  know  how  sound  propagates  in  the  ocean  as  a 
function  of  frequency  for  different  source/receiver  configurations 
and  for  different  environmental  conditions. 

By  now  the  theory  of  acoustic  propagation  is  well  developed, 
providing  both  a  good  general  understanding  and  a  detailed 
description  of  how  sound  travels  in  the  ocean.  The  theoretical 
basis  is  the  acoustic  wave  equation,  which  has  to  be  solved  with 
realistic  boundary  conditions  at  the  sea  surface  and  at  the  sea 
floor.  This  problem  is  generally  too  complex  for  applying  analy- 
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tical  solutions,  and  hence  we  must  resort  to  numerical  methods. 
Several  solution  techniques  (ray,  mode,  FFP,  PE)  have  been  intro¬ 
duced  over  the  years,  with  the  acoustic  models  increasing  In 
complexity  as  computers  became  faster  and  more  powerful.  Ray 
theory  was  the  only  practical  technique  for  solving  propagation 
problems  until  the  beginning  of  the  1970s.  Then  advances  in  com¬ 
puter  technology  made  it  possible  to  consider  solving  more  complex 
equations,  and,  consequently,  new  techniques  (normal  mode,  fast 
field,  parabolic  equation)  came  into  extensive  use  during  the  last 
decade. 

■*-  In  this  paper  we  outline  the  basics  of  sound  propagation  in 
the  ocean,  including  important  propagation  and  loss  mechanisms, 
and  a  simplified  environmental  description  for  use  in  numerical 
models.  We  then  proceed  to  outline  the  mathematical  foundation  of 
the  most  widely  used  numerical  models,  and  we  demonstrate  the  con¬ 
sistency  and  inter-relationship  between  the  various  models  through 
a  few  numerical  examples  for  deep  and  shallow  water  environments. 
We  then  indicate  the  areas  of  applicability  of  the  various 

acoustic  models  taking  into  account  both  limitations  in  the 

underlying  theory  and  the  numerical  efficiency  of  the  actual  com¬ 
puter  codes.  A  demonstration  of  the  wide  range  of  applicability 
of  these  ocean-acoustic  models  is  provided  by  a  series  of  com¬ 
putational  examples,  where  the  models  have  been  applied  to  some 

general  wave  propagation  problems,  including  beam  reflection  at 
fluid/solid  Interfaces,  propagation  from  ducts  into  free-space, 
up-  and  down-slope  propagation  involving  mode  coupling  and  mode 
cutoff. 


2  SOUND  PROPAGATION  IN  THE  OCEAN 

The  goal  of  ocean  acoustic  modelling  Is  to  estimate  the  spa¬ 
tial  properties  of  the  sound  pressure  'field  as  a  function  of 
source  frequency.  To  clarify  the  complexity  of  the  modelling 
problem,  let  us  briefly  review  the  environmental  acoustics  of  the 
ocean.  Figure  1  is  a  schematic  of  some  important  sound  propaga¬ 
tion  paths;  two  possible  sound-source  locations  are  cr.  the  left 
and  sound  is  propagating  from  left  to  right.  Two  dashed  lines  at 
0  and  80  km  range  indicate  two  of  the  innumerable  ways  in  which 
sound  speed  in  the  water  can  vary  with  depth  from  place  to  place 
(or  from  time  to  time).  Lines  A,  B,  C,  and  D  represent  four 
possible  sound-propagation  paths  whose  shapes  are  determined  by 
the  location  of  the  source  and  the  sound-speed  structure  over  the 
extent  of  the  propagation. 

Path  A  from  the  shallow  source  is  "surf are-duct ”  propagation, 
because  the  sound-speed  profile  is  such  that  the  sound  is  trapped 
near  the  surface  of  the  ocean.  Paths  B,  C,  and  D  are  from  the 
deeper  source.  Ray  B,  leaving  the  source  at  a  small  angle  from 
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RANGE  (km) 


Fig.  1  Schematic  of  sound  propagation  in  the  ocean 


the  horizontal,  tends  to  propagate  in  the  "deep  sound  channel" 
without  interacting  with  the  boundaries  (surface  and  bottom)  of 
the  ocean.  At  slightly  steeper  angles  (path  C)  we  have 
"convergence  zone"  propagation,  which  is  a  spatially  periodic  phe¬ 
nomenon  of  zones  of  high  intensity  near  the  surface.  Here  sound 
interacts  with  the  ocean  surface  but  not  with  the  bottom.  Path  D 
is  the  "bottom-bounce  path",  which  has  a  shorter  cycle  period  than 
the  convergence  zone  path.  The  right-hand  side  of  Fig.  1  depicts 
propagation  on  the  continental  shelf  (shallow  water)  where  a 
complicated  bottom  structure  combined  with  variable  sound-speed 
profiles  result  in  rather  complicated  propagation  conditions  not 
always  suited  for  a  ray  representation. 

Our  ability  to  model  acoustic  propagation  effectively  in  the 
ocean  is  determined  by  the  accuracy  with  which  acoustic  loss 
mechanisms  in  the  ocean  environment  are  handled.  Aside  from 
geometrical  spreading  loss  (spherical,  cylindrical,  etc.)  the  main 
loss  mechanisms  are  volume  absorption,  bottom-reflection  loss, 
surface  and  bottom  scattering  loss. 

Volume  absorption  in  sea  water,  caused  by  viscosity  and  che¬ 
mical  relaxation,  increases  with  increasing  frequency.  This  loss 
mechanism  is  the  dominant  attenuation  factor  associated  with  path 
B  in  Fig.  1,  since  this  path  does  npt  interact  with  the  boun¬ 
daries.  Because  there  is  very  little  volume  absorption  at  low 

frequencies,  deep-sound-channel  propagation  has  been  observed  to 
distances  of  many  thousands  of  kilometres. 

When  sound  Interacts  with  the  sea  floor,  the  nature  of  the 

bottom  becomes  important.  Figure  2  depicts  simple  bottom-loss 
curves,  with  zero  dB  loss  indicating  perfect  reflection.  For  an 
ideal  bottom  without  volume  absorption  (non-lossy)  we  still  get 
severe  reflection  loss  above  a  certain  critical  angle  0C  due  to 

transmission  into  the  bottom.  For  a  real  bottom  with  volume 
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GRAZING  ANGLE  [deg) 

Fig.  2  Bottom  loss  versus  grazing  angle 

absorption  (lossy)  we  never  get  perfect  reflection,  even  though 
the  curves  look  similar.  Path  D  In  Fig.  1,  the  bottom-bounce 
path,  often  corresponds  to  angles  near  or  above  the  critical 
angle;  therefore  after  a  few  bounces  It  is  highly  attenuated.  On 
the  other  hand,  for  shallower  angles,  many  more  bounces  are 
possible;  hence  in  shallow  water  (path  E)  most  of  the  energy  that 
propagates  is  close  to  the  horizontal.  In  reality,  much  of  the 
ocean  bottom  Is  layered  and  also  supports  shear  waves;  in  this 
case  bottom  loss  becomes  a  complicated  function  of  frequency  and 
grazing  angle.  The  overall  effect  of  bottom  loss  on  sound  propa¬ 
gation  in  the  water  column  is  an  increasing  loss  with  decreasing 
frequency. 


So 


Fig.  3  Environmental  input  to  ocean  acoustic  models 
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A  rough  sea  surface  or  sea  floor  causes  scattering  of  the 
Incident  sound.  The  result  Is  a  decay  of  the  mean  acoustic  field 
in  the  water  column  as  a  function  of  range  (scattering  loss),  with 
the  scattered  energy  being  lost  to  the  ocean  bottom  through  steep- 
angle  propagation.  The  scattering  loss  increases  with  increasing 
frequency,  and  the  propagation  paths  mainly  affected  are  paths  A 
and  C  (surface  scattering  loss)  and  paths  D  and  E  (surface  and 
bottom  scattering  loss). 

A  consistent  mathematical  model  of  sound  propagation  in  the 
ocean  must  contain  the  physics  that  governs  the  above-mentioned 
propagation  and  loss  mechanisms.  A  summary  of  the  environmental 
Inputs  needed  for  a  realistic  description  of  the  ocean  waveguide 
is  given  in  Fig.  3.  In  this  simplified  model  the  ocean  consists 
of  a  water  column  of  depth  H0  limited  by  a  rough  sea  surface  and  a 
rough  sea  floor.  The  sound  speed  c0  in  the  water  column  may  vary 
arbitrarily  with  depth,  while  density  p0  and  attenuation  0O  are 
considered  constant.  Even  though  real  ocean  bottoms  exhibit  a 
complicated  layering,  we  have  found  that  a  simple  two-laye1"*'1 
geoacoustic  model  generally  provides  the  necessary  degr 
freedom  to  accurately  include  bottom  effects  in  numerica  .els 
for  many  ocean  areas.  Hence  the  bottom  may  consist  ol  just  a 
sediment  layer  of  thickness  Hj  and  a  semi-infinite  subbottom.  T  ' 
model  should  allow  for  sound  speed,  density,  and  attenuatiu  j 
vary  arbitrarily  with  depth  in  the  sediment  layer,  while  the  sub¬ 
bottom  can  be  considered  homogeneous.  It  is  desirable  that  the 
model  can  handle  shear-wave  propagation  in  both  bottom  layers. 
Finally,  in  real  ocean  environments  the  parameters  given  in  Fig.  3 
may  all  vary  with  range. 

3  MATHEMATICAL  FOUNDATION  OF  OCEAN . ACOUSTIC  MODELS 

We  briefly  present  the  mathematical  foundations  of  the  four 
models  discussed  in  this  paper.  A  more  detailed  description  can 
be  found  in  references  <1-10>. 

The  starting  point  for  all  the  models  is  the  wave  equation 
for  a  harmonic  point  source  with  time  dependence  exp(-iwt), 

v2$(x,y,z)  +r~7~ - r]  4>(x,y,z)  =  -6(x-x0)6(y-y0)6(z-z0) 

Lc(x,y,z)J  (1) 

*  $exp(-i«t)  (2) 


At  any  point 
satisfies  Eq.  (1) 


(x,y,z)  in  the  medium,  the  velocity  potential  <J> 
where  c(x,y,z)  is  the  sound  speed  of  the  medium 
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and  6  is  the  Dirac  delta  function.  The  source  is  at  the  coor¬ 
dinate  (x0,y0,z0)  where  z  is  the  depth  coordinate,  which  is  taken 
to  be  positive  in  the  downward  direction  from  the  ocean  surface. 


For  the  boundary  condition  at  the  ocean  surface  we  take  the 
density  of  air  to  be  negligible  compared  with  that  of  water; 
hence,  the  pressure  must  vanish  at  the  ocean  surface 
( "pressure-release  surface").  At  a  boundary  between  two  media 
such  as  the  ocean  and  the  ocean  bottom,  the  balancing  of  forces  at 
the  Interface  require  that  physical  quantities  such  as  particle 
velocity  and  pressure  be  continuous  across  the  boundary: 


x,y,  or  z 


(3) 


p  “  -i(vp$ 


(4) 


If  the  ocean  bottom  is  treated  as  an  elastic  medium  that  can 
support  shear  motions,  there  is  the  additional  boundary  condition 
that  tangential  stress  must  be  continuous.  Since  the  water  column 
cannot  support  shear  waves,  this  requires  that  the  tangential 
stress  In  the  ocean  bottom  vanishes  at  the  interface. 

Four  widely  used  solution  techniques  for  Eq.  (1)  are  schema¬ 
tically  represented  in  Fig.  4.  The  derivation  of  the  classical 
ray  solution  can  be  found  in  most  text  books  on  acoustics,  as  can 
the  details  of  the  well-established  normal-mode  solution.  The 
fast-field  technique  is  not  yet  in  standard  text  books,  but  it  is 
a  powerful  tool  for  solving  propagation  problems  in  stratified 
environments.  The  parabolic  equation  technique  is  a  recent  advent 
in  acoustic  modelling.  This  method  is  particularly  suited  for 
propagation  in  range-dependent  environment's. 

* 

We  briefly  describe  the  derivation  of  the  above  four  solution 
techniques,  starting  with  range-independent  wa  'e  theory. 


WAVE  EQUAT ION 


Fig.  4  Four  techniques  for  solving  the  wav'*  egra*"'on 
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(39) 


where  l  is  the  arc  length  along  a  ray  and  X  is  the  coordinate.  It 
can  be  shown  that  the  direction  of  the  average  energy  flow  is 
along  these  trajectories  and  the  amplitude  of  the  field  at  any 
point  can  be  obtained  from  the  density  of  these  rays;  formally, 
having  solved  for  S,  the  amplitude  is  obtained  from  solving  the 
second  part  of  Eq.  (36).  We  also  mention  here  that  corrected  ray 
theory  assumes  that  iji  is  a  function  of  frequency  and  an  expansion 
in  powers  of  inverse  frequency  is  made,  the  leading  term  being  the 
in£jL  nite-f  requency  solution  with  the  additional  terms  being 
corrections  from  the  infinite-frequency  solution. 


The  advantages  of  ray  theory  methods  are  that  the  com¬ 
putations  are  rapid  and  that  ray  traces  give  a  very  physical  pic¬ 
ture  of  the  acoustic  paths.  The  disadvantage  is  that  ray-theory 
is  an  infinite-frequency  approximation  and  therefore  does  not 
include  diffraction  and  other  wave  effects.  This  shortcoming  also 
prevents  ray  theory  from  adequately  describing  significant  bottom 
interaction  and  low-frequency  ducted  propagation. 


4  NUMERICAL  MODELS;  THEIR  APPLICABILITY  AND  CONSISTENCY 

The  four  acoustic  models  described  in  this  paper  are  a  repre¬ 
sentative  subset  of  the  many  different  propagation  models  in  use 
in  underwater  acoustics  today.  The  reason  for  developing  new 
models  is  to  obtain  either  more  accurate  solutions  or  faster  solu¬ 
tions  to  specific  problems.  Each  model  has  its  area  of  applicabi¬ 
lity  depending  on  the  theoretical  limitations  in  the  model  and  on 
the  numerical  efficiency  of  the  computer  code. 


MODEL  TYPE 

APPLICATIONS 

SHALLOW  WATER 

DEEP  WATER 

LF _ 

_ hf 

LF 

HF 

wm 

RD 

RI 

RD 

RI 

RD 

RI 

RAY 

Hi 

NORMAL  MODE 

Hi 

|p 

Wk 

Wk 

'm 

FAST  FIELD  (FFP) 

PARABOLIC  EQ.(PE) 

'Mb 

I8S 

m 

LF:  LOW  FREQUENCY  (<  500  HI)  RJ:  RANGE  - 1 NOEPENDEN T  ENVIRONMENT 

HF:  HIGH  FREQUENCY  <>  500  HI)  R D :  RANGE-DEPENQENT  ENVIRONMENT 

Fig.  9  Applicability  of  four  propagation  models 
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The  parabolic  wave  equation  as  given  by  Eq.  (29)  is  based  on 
the  paraxial  approximation,  and  hence  only  propagation  close  to 
the  horizontal  (±  20°)  is  accurately  handled.  This  angle  limita¬ 
tion  is  of  minor  importance  for  a  wide  class  of  ocean  acoustics 
problems.  However  in  studying  bottom-interacting  propagation,  the 
narrow-angle  approximation  becomes  a  serious  limitation.  It  has 
recently  been  shown  <35,36>  that  a  slight  modification  to  Eq.  (20) 
can  improve  the  angle  coverage  to  ±  40°,  yielding  a  modified  para¬ 
bolic  equation,  which,  however,  can  no  longer  be  solved  by  the 
split-step  technique.  Instea'1  finite-difference  solution  tech¬ 
niques  have  been  applied  <35,37>,  and  a  working  computer  code  is 
already  available  <38>. 

Numerical  PE  results  given  in  this  paper  were  all  done  with  a 
model  based  on  the  standard  parabolic  equation  technique  as  deli¬ 
neated  above. 


3.4  Ray  theory 


This  paper  is  concerned  mainly  with  wave  theory;  neverthe¬ 
less,  for  completeness,  we  include  a  brief  description  of  ray 
theory.  In  this  case  we  assume  a  solution  of  Eq.  (1)  (with  right- 
hand  side  equal  to  zero)  as 


♦  -  <Kx,y,z) 


iS(x,y,z) 

6  t 


(35) 


S(x,y,z)  is  a  phase  function  that  includes  rapid  variations  as  a 
function  of  range,  and  ♦(x,y,z)  is  a  more  slowly  varying  envelope 
function  in  which  geometrical  spreading  and  loss  mechanisms  are 
Included  (in  the  PE,  S  contains  the  cylindrical  spreading  factor). 
Substituting  Eq.  (35)  into  the  wave  equation  and  separating  real 
and  imaginary  parts,  we  obtain 


i-V2<|<  -  (VS)2  +  k2  =  0,  2( Vi|/  •’  VS)  +  t|/V2S  =  0. 

We  now  make  the  geometrical-acoustics  approximation 

£  V2*  «  k2  , 


(36) 


(37) 


that  is,  the  amplitude  of  the  phase  function  varies  slowly  in 
range  with  respect  to  wavelength.  Substituting  Eq.  (37)  in  Eq. 
(36)  gives  the  eikonal  equation, 

(VS)2  -  k2.  (38) 

The  trajectory  of  the  rays  is  perpendicular  to  the  surfaces  of 
constant  phase  (wavefronts),  S,  and  is  expressed  by 
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RANGE  ( KM  ) 


Fig.  8  Computed  transmission  losses  for  test  problem 


In  Sect.  3.2  we  compared  FFP  and  NM,  using  a  simple  environ¬ 
ment  (Fig.  5)  as  an  example.  We  now  look  at  transmission  loss 
from  the  point  of  view  of  discrete,  dlscrete-plus-contlnuous,  and 
PE,  which  does  not  obviously  distinguish  between  regions  in  the 
spectrum.  Figure  8  shows  computed  transmission  loss  from  0  to  6 
km,  defined  as  TL  *  -20  log  (P/Pi),  wfrere  Pj  Is  the  pressure 
amplitude  at  1  m  distance  from  the  source.  Note  that  the  PE 
tracks  the  FFP  results  In  the  nearfield  indicating  that  at  least 
part  of  the  continuous  portion  of  the  spectrum  Is  Included  In  the 
PE  calculation  when  using  a  gaussian  starting  field.  We  can  also 
see  how  the  three  model  results  converge  In  the  farfield;  recall 
that  the  NM  calculation  does  not  Include  the  full  nearfield 
contribution.  This  particular  example  clearly  illustrates  the 
consistency  and  Inter-relationship  between  the  three  models. 

The  advantages  of  the  PE  are  that  It  handles  a  range- 
dependent  environment  and  gives  the  acoustic  field  in  the  entire 
water  column  without  additional  computational  effort.  Its  disad¬ 
vantages  are  that  the  procedure  Is  not  easily  automated,  and  it  is 
practical  only  for  low-frequency  propagation  since  computation 
times  Increase  with  frequency  squared.  Moreover,  there  is  no 
straightforward  way  of  handling  shear  propagation  in  the  bottom. 
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m 

>l>(r,z)  -  jf  iKr0,s)  •  e 


where  Ar  ■*  r  -  r0. 

By  introducing,  the  symbol  s  for  the  fourier  transform  from  the  z- 
domain  and  as  the  inverse  transform,  Eq.  (33)  may  be  written 

as 

Hi®.  (n2-l)Ar  _  s2 

♦  <r  +  Ar.z)  -  e  2  •  S  *  /  e  2k°  '  ,z) }) . 

(34) 

Equation  (34)  is  the  so-called  "split-step”  marching  solution 
of  the  parabolic  equation.  The  fourier  transforms  are  performed 
using  an  FFT.  It  is  the  solution  for  n  constant,  but  the  error 
Introduced  when  n  (profile  and  bathymetry)  varies  with  range  and 
depth  can  be  made  arbitrarily  small  by  increasing  the  transform 
size  and  decreasing  the  range-step  size  <33,34>. 

The  parabolic  equation  is  not  a  boundary-value  equation  as  we 
have  numerically  formulated  it  above.  We  can  Include  the  surface 
boundary  condition  by  taking  an  anti-symmetric  FFT  about  the  sea 
surface  (z  ■  0).  In  practice  this  is  performed  by  taking  sine 
transforms.  The  boundary  conditions  in  the  bottom  are  simulated 
by  including  the  discontinuity  in  velocity  in  the  sound-speed  pro¬ 
file.  There  are  methods  to  also  include  the  density  discontinuity 
<33>.  The  radiation  condition  as  z  goes  to  Infinity  is  simulated 
by  requiring  the  field  to  exponentially  tail  off  for  large  values 
of  z  beyond  which  there  would  not.  .be  any  significant  acoustic 
interaction. 

As  mentioned  above,  the  PE  method  requires  an  initial 
starting  solution.  Two  methods  have  been  used  for  describing  a 
point  source.  The  first  method  is  to  initialize  the  field  with  a 
set  of  normal  modes  descriptive  of  the  point  source  in  the 
starting  environment.  This  would  not  include  the  continuous  por¬ 
tion  of  the  spectrum  (see  Sect.  3.2),  but  for  long-distance  propa¬ 
gation  this  approximation  is  adequate.'  A  second  approach  has 
proved  to  be  simpler  and  as  effective.  The  point  source  is 
approximated  by  two  gaussians  that  are  anti-symmetric  about  the 
sea  surface,  thereby  automatically  including  the  pressure-release 
boundary  condition  at  the  surface.  Both  starting  techniques  have 
been  used  in  this  paper.  We  will  see  that  by  using  the  gaussian 
starting  field  part  of  the  continuous  spectrum  is  included  in  the 
PE  solution. 
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This  particular  example  shows  the  consistency  between  FFP  and 
NM  calculations.  In  the  next  section  we  show  the  results  on 
transmission  loss  when  discrete  and  dtscrete-plus-contlnuous 
spectra  are  Included  and  we  compare  them  with  calculations  using 
the  parabolic  equation  technique. 


3.3  Parabolic  equation 

If  the  environment  varies  both  in  range  and  depth,  the  wave 
equation  cannot  be  separated  and  therefore  direct  numerical 
integration  is  required.  At  present  there  are  no  practical 
methods  to  perform  this  direct  integration  of  the  three- 
dimensional  wave  equation,  which  is  a  boundary-value  problem.  An 
alternative  approach  is  to  derive  an  approximate  wave  equation 
that  lends  Itself  to  practical  numerical  solution.  We  now  outline 
the  derivation  of  such  an  approximation,  the  Parabolic  Equation 
(PE),  which  was  introduced  into  underwater  acoustics  in  1973  by 
Tappert  and  Hardin  <17,33>. 


The  velocity  potential  is  decomposed  as  follows: 

♦  -  <|>(r ,z)  •  S(r),  (23) 

and  we  substitute  $  into  Eq.  (1)  in  a  source-free  region: 


3!i+ii!  + 

3r2  3z2 


Liii 

S  3rJ 


“0. 

(24) 


That  is,  we  will  eventually  end  up  with  an  equation  that 
allows  a  "marching-in-range ”  solutlpn  and  we  will  have  to  ini¬ 
tialize  the  solution  in  some  way  (see  below).  We  use  the  notation 
that 


k2  - 


k2n2 , 


(25) 


where  n  is  an  "index  of  refraction”  equal  to  cG/c,  where  cQ  is  a 
reference  speed. 

Equation  (24)  may  be  separated  into  two  differential 
equations  by  setting  the  terms  in  the  first  bracket  equal  to  -Sk2 
and  the  terms  in  the  second  bracket  equal  to  <|<k2  ,  where  k2  is 
the  separation  constant.  The  functions  S(r)  and  <|>(r,z)  then  have 
to  satisfy  the  following  two  equations: 


32S 

3r2 


+ 


r  3r  o 


0 


(26) 
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TABLE  1 

MODAL  EIGENVALUES 


Mode  no. 

Wavenumber  (m~l) 

1 

0.10423 

2 

0.10386 

3 

0.10356 

4 

0.10328 

5 

0.10297 

6 

0.10261 

7 

0.10221 

8 

0.10182 

9 

0.10142 

10 

0.10102 

11 

0.10064 

Figure  7  displays  the  amplitudes  of  the  eleven  inodes  plotted 
as  a  function  of  depth.  The  dashed  line  indicates  the 
source/receiver  position.  Notice  the  high  excitation  of  the 
second  mode  and  the  low  excitation  of  the  third,  sixth,  and  ninth 
modes,  and  notice  the  one-to-one  correspondence  with  the  FFP 
integrand  shown  in  Fig.  6b. 


Fig.  7  Mode-amplitude  functions  for  test  problem 
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Fig.  6  FFP  integrand  for  test  problem 

a)  full  spectrum,  b)  discrete  spectrum 
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We  shall  demonstrate  the  inter-relationship  between  normal 
mode  (NM)  and  fast  field  (FFP)  solutions  through  a  numerical 
example  using  the  environment  given  In  Fig.  5.  The  upward- 
refracting  sound-speed  profile  defines  a  surface  duct  of  thickness 
1500  m.  We  consider  propagation  for  source  and  receiver  both  a 
500  m  depth,  and  for  a  frequency  of  25  Hz.  Two  plots  of  the  FFP 
integrand  are  shown  in  Fig.  6.  The  horizontal  axis  Is  the  wave- 
number  (n)  and  the  vertical  axis  is  a  normalized  absolute  value  of 
the  integrand  amplitude.  The  vertical  dashed  line  in  Fig.  6a 
separates  the  discrete  spectrum  (defined  by  Ineq.  18)  and  the  con¬ 
tinuous  spectrum.  A  blowup  of  the  discrete  spectrum  is  shown  in 
Fig.  6b.  In  this  particular  case  there  are  eleven  propagating 
mo<tes.  The  asterisks  on  the  plot  are  the  locations  and  the  ampli¬ 
tudes  of  the  modes  from  a  NM  calculation  <18>.  We  see  that  the 
eigenvalues  as  calculated  from  the  NM  model  (Table  1)  coincide 
with  the  peaks  (poles)  in  the  FFP  integrand,  and  also  that  the  NM 
amplitudes  correspond  to  the  amplitudes  of  the  peaks.  (It  turns 
out  that  this  one-to-one  correspondence  in  the  amplitudes  is 
because  there  is  virtually  no  loss  in  this  problem.  Otherwise, 
the  correspondence  would  not  be  as  precise  because  loss  shows  up 
in  the  FFP  calculation  as  widths  in  the  peaks;  nevertheless,  for 
realistic  losses  the  location  of  the  poles  would  be  the  same). 


Fig.  5  Sound-speed  profile  for  test  problem 
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ip(z0)  -i»/4  un(z0)un(z)  iknr 
-  e  E  -  e 

(8wr)V2  n  k^2 


(22) 


In  addition  to  the  decay  of  the  field  due  to  cylindrical 

spreading,  other  loss  mechanisms  such  as  volume  attenuation  in  the 
water  column  and  bottom  are  Included  in  Eq.  (22)  because  the 
eigenvalues,  kj,,  have  positive  imaginary  parts  <16>,  thereby 
resulting  in  an  exponential  attenuation  of  each  normal-mode  term. 
Equation  (22)  gives  us  the  important  result  that  the  field  at  a 

depth  z  is  proportional  to  a  sum  of  the  products  of  normal  modes 
evaluated  at  the  source  and  the  receiver  depth.  The  normal  modes 
are  the  "natural  vibrations"  of  the  system  and  if  a  point  source 
is  located  at  the  null  of  a  particular  normal  mode,  that  mode  will 
not  be  excited.  Similarly,  if  a  point  receiver  is  placed  at  the 
null  of  a  particular  mode,  that  mode  contribution  to  the  total 
field  will  not  be  sensed. 

In  analogy  to  the  FFP  procedure,  the  main  numerical  effort 
for  the  normal-mode  procedure  is  the  solution  of  the  eigenvalue 

problem  defined  by  Eq.  (16)  and  the  boundary  conditions.  There 
are  many  techniques  to  solve  this  equation  <17  — 18>  but  they  are 
mainly  applicable  to  low-frequency  or  shallow-water  propagation, 
where  there  is  only  a  small  number  of  modes  <19>.  However,  there 

are  also  techniques  to  handle  deep-water  high-frequency  propaga¬ 
tion  using  normal  modes  <20>. 

The  advantages  of  the  normal-mode  procedure  are,  first,  that 
once  Eq.  (16)  is  solved  we  have  the  solution  for  all  source/ 
receiver  configurations,  and,  second,  that  the  whole  solution  pro¬ 
cedure  can  be  highly  automated.  In  addition,  the  normal-mode  pro¬ 
cedure  can  be  easily  extended  to  slightly  range-dependent 
environments  using  the  adiabatic  approximation  where  mode  coupling 
is  neglected  <21— 26> ;  numerical  methods  for  Including  mode- 
coupling  effects  are  present  areas  of  research  <27—3 1> .  The 
disadvantages  of  the  normal-mode  solution  are  that  conventional 
procedures  do  not  include  nearfleld  effects  (the  exception  is 
Stickler's  work  <32>,  but  even  there  the  nearfield  is  evaluated 

with  a  procedure  similar  to  the  FFP  approach)  and  there  are 
restrictions  on  how  one  can  treat  shear  propagation  in  the 

bottom. 

Both  FFP  and  normal  modes  are  solutions  of  Eqs.  (5)  and  (6). 
The  difference  between  the  two  is  that  the  normal-mode  method 
restricts  the  integration  to  horizontal  wavenumbers  in  the  inter¬ 
val  corresponding  to  the  discrete  portion  of  the  spectrum  defined 
by  Ineq.  (18).  From  Eq.  (20),  we  see  that  the  integrand  has 

poles  at  n  ■  kjj.  Hence,  the  function  u  in  Eq.  (10)  and  Xm  of  Eqs. 

(13)  and  (14)  should  have  poles  at  the  same  locations. 
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tion,  we  require  the  un(z)  be  bounded  as  z  ♦  “.  The  normal  modes 
un(z)  form  a  complete  orthogonal  set  that  satisfies  the  relation 


/  p(z)un(z)um(z)dz  -  inm, 
o 


07) 


where  the  density  p(z)  takes  its  appropriate  value  in  each  layer 
and  6nm  is  the  Kronecker-delta  symbol.  The  spectrum  of  eigen¬ 
values  consists  of  a  discrete  part  and  a  continuous  part,  the 
discrete  eigenvalues  occurring  In  the  interval 

*■  w/c2  <  k,,  <  maxf<o/c(z)  J ,  (18) 

where  C2  is  the  highest  speed  of  the  system.  In  the  present 
treatment  we  consider  only  the  discrete  eigenvalues,  since,  in 
general,  the  continuous  spectrum  makes  a  negligible  contribution 
beyond  the  nearfield  of  the  source  (and  requires  an  FFP-type 
calculation  in  any  event). 


We  now  substitute  Eq.  (15)  into  Eq.  (6),  multiply  the 
resulting  equation  by  p(z)um(z),  and  integrate  over  z  from  0  to 
giving: 


an 


1  P(z0)un(z0) 


n 


(19) 


Substituting  Eqs.  (15)  and  (19)  back  into  Eq.  (5)  we  obtain  an 
Integral  representation  of  the  velocity  potential, 


p(z  )  m  00 

♦  (x.y.z)  “  — 2-  /  dnx  /  dny  z 

(2*)2  -•  n 

*  exp[i(nxx  +  nyy)J. 


un(z0)»n(z) 

*  •  n2-kz 
n 


We  evaluate  the  above  integral  by  choosing  a  path  about  the 
poles  so  as  to  lead  to  an  outgoing  wave  from  the  source  point 
r  -  0.  Each  Integral  in  Eq.  (20)  is  proportional  to  the  two- 
dimensional  plane-wave  representation  of  the  zero-order  Hankel 
function  of  the  first  kind  <15>,  and  therefore  $(x,y,z)  can  be 
expressed  as: 

♦  (r,z)  -  ~p(z0)  l  un(z0)un(z)Hol^(knr).  (21) 

4  n 


The  asymptotic  form  of  the  Hankel  function  can  then  be  used  to 
obtain 
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putatlon.  Numerical  procedures  for  Integrating  Eq.  (6)  are  given 
in  references  <11-13>. 

The  first  implementation  of  the  FFP  approach  was  done  by 
DiNapoll  <11>  around  1970.  His  method  is  for  a  stratified  f luid 
environment,  and  the  numerical  procedure  is  quite  fast,  since  it 
uses  recurrence  relations  of  hypergeometric  functions  to  solve  Eq. 
(6)  for  all  the  wavenumbers  (n ’s)  rather  than  solve  the  equation 
for  one  n  at  the  time.  A  more  general  solution  technique  was 
devised  by  Kutschale  <12>  for  an  arbitrary  stratification  of  solid 
layers.  This  technique  employs  the  Thoms on-Haske 11  matrix  method, 
and  Eq.  (6)  is  here  solved  separately  for  each  wavenumber.  This 
solution  technique  is  computationally  quite  slow  and  there  is  no 
efficient  way  of  doing  calculations  simultaneously  for  many  sour¬ 
ces  and  receivers.  The  most  recent  FFP  technique  was  developed  by 
Schmidt  <13>.  Again  the  solution  is  for  an  arbitrary  stratifica¬ 
tion  of  homogeneous  solid  layers.  Displacements  and  stresses  are 
expressed  in  terms  of  three  scalar  potentials  for  each  layer,  as 
described  in  <14>.  Boundary  conditions  are  then  matched  at  each 
interface  yielding  a  linear  system  of  equations  in  the  Hankel 
transforms  of  the  potentials.  Equation  (6)  is  again  solved  at 
discrete  horizontal  wavenumbers;  with  the  coefficient  matrices 
being  of  band  form,  the  equations  are  solved  very  efficiently  by 
gaussian  elimination.  This  solution  technique  is  considerably 
faster  than  the  Kutschale  technique,  and  it  furthermore  allows  for 
an  efficient  evaluation  of  the  acoustic  field  for  many  sources  and 
receivers  at  a  time.  This,  in  turn,  means  that  the  model  can  be 
applied  to  a  variety  of  new  problems,  including  beam  reflection 
problems  as  shown  in  Sect.  5.1. 

The  main  advantage  of  the  FFP  is  that  it  provides  the  full 
solution  to  sound  propagation  in  a  multilayered  solid  medium,  and 
hence  consitutes  a  benchmark  against  which  other  approximate  solu¬ 
tions  can  be  checked.  Its  main  disadvantage  is  that  the  procedure 
is  not  easily  automated. 

3.2  Normal-mode  solution 


The  alternative  to  a  direct  numerical  integration  of  Eq.  (6) 
is  to  expand  u  into  a  complete  set  of  normal  modes: 

■  E  an^r|x*rly)un(z)  ,  (15) 

where  the  t^'s  are  the  solutions  to  the  eigenvalue  equation 

- ---—  +  fk2(z)  -  k2]un(z)  -  0  (16) 

dz2 

that  satisfies  the  above-mentioned  boundary  conditions.  In  addi- 
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greater  than  a  few  wavelengths  from  the  source,  the  asymptotic 
form  of  the  Hankel  function  can  be  used: 

H^(nr)  ~  (2/wnr)1^2  exp[  'nr  -  ir/4)] 


and,  hence,  Eq.  (9)  can  be  expressed  as 


4>(r,z) 


e-iir/4 


~  J  dn  /nu(n,z)elnr. 


(10) 


where  the  factor  1//F  indicates  cylindrical  spreading. 

Equation  (10)  can  be  numerically  integrated  to  obtain  the 
acoustic  field  at  the  range  r  and  depth  z.  In  order  to  do  this  we 
must  solve  Eq.  (6)  for  many  n's  to  have  a  sufficient  set  of  u's 
as  a  function  of  n  so  that  the  integration  over  n  in  Eq.  (10)  can 
be  performed.  Given  that  u  has  been  obtained  numerically  as  a 
function  of  n,  the  integration  can  be  done  using  an  FFT  algorithm. 
This  total  procedure  is  called  the  Fast  Field  Program  (FFP) 
<11-13>,  although  most  of  the  numerical  effort  goes  into  solving 
Eq.  (6)  for  the  many  n's.  For  computation  we  discretize  Eq.  (10) 
by  letting 


nm  ”  nc  +  m&n;  rn  -  r0  +  nAr;  (m,n)  =  0,1 ,2... ,N-1 , 

with  the  additional  relation 

ArAn  =  — 

N 


(11) 


(12) 


and  N  is  an  integral  power  of  two.  Note^  that  the  discretization 
relations  of  Eq.  (11)  restricts  the  solution  to  outgoing  waves. 
Substituting  Eq.  (11)  into  Eq.  (10)  we  obtain 


4>(rn,z) 


An 


lOo^n 

e 


N-l 

l 

W0 


Xm 


i2vmn/N 

e 


(13) 


and  hence  the  input  to  the  FFT  is 

. —  imr,jAn 

Xm  “  *'nm  u(nm,z)e 


(14) 


Equations  (13)  and  (14)  specify  the  numerical  procedure  to  be 
employed  in  solving  the  wave  equation  using  the  FFP  approach  after 
Eq.  (6)  has  been  solved  numerically  for  the  complete  set  of  u's  as 
a  function  of  n.  As  mentioned  above,  the  main  effort  in  the  FFP 
approach  is  the  numerical  integration  of  Eq.  (6)  and  not  the  final 
implementation  of  Eqs.  (13)  and  (14),  which  is  a  simple  FFT  com- 
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3.1  Fast  field  solution 

Here  we  are  solving  the  wave  equation  for  the  case  where  the 
sound-speed  profile  Is  only  a  function  of  depth  and  the  bottom  Is 
flat;  this  type  of  environment  Is  often  referred  to  as  the  hori¬ 
zontally  stratified  ocean.  From  Eq.  (1)  we  therefore  have  that 
c(x,y,z)  Is  simply  c(z).  Because  the  environment  Is  Independent 
of  “r",  the  horizontal  coordinates  (x,y),  one  possible  method  of 
solving  Eq.  (1)  Is  to  Fourier  decompose  the  acoustic  field  into  an 
Infinite  set  of  horizontal  waves: 

♦  (x,y,z)  -  /  d2n  u(n,z)  eln*r  .  (5) 

Substituting  Eq*  (5)  Into  Eq.  (1)  we  obtain  the  equation  for 
t  z)  , 

- — 2 - +  lk2(z)  -  n2]u(nXfny;z)  -  -  ^-6(z-z0),  (6) 

dz4 

where  k(z)  -  u/c(z)  and  n2  -  n2  +  n2  ,  is  the  horizontal  wave- 
number  of  the  individual  plane  waves. 

Using  polar  coordinates  we  can  rewrite  Eq.  (S)  as 

♦ (r , z)  -  A-  J*  d«  /"  n  dn  u(n,z)elnrcose.  (7) 

2xo  o 

We  now  integrate  over  the  azimuthal  angle  to  obtain 
00 

♦(r,z)  -  /  n  dn  u(n,z)J0(nr),  «  (8) 

o 

where  JQ  is  the  zero1*1  order  Bessel  function.  Using  the  rela¬ 
tionship  that 

J0(nr)  -  1/2  [H<°(nr)  +  H^^nr)], 

where  the  H's  are  Hankel  functions  and  noting  from  Eq.  (6)  that 
u(n,z)  is  even  in  n,  we  can  rewrite  Eq.  (8)  as 

♦  (r,z)  -1/2  /  U  dn  u(n,z)  H^(nr),  (9) 

—00  ® 

where  now  the  integration  over  n  is  from  -®  to  «.  For  ranges 
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To  Indicate  with  some  precision  the  type  of  ocean  environment 
for  which  a  Riven  model  should  be  used,  we  have  classified 
environments  according  to  water  depth,  frequency,  and  environmen¬ 
tal  complexity,  as  shown  in  Fig.  9.  Here  shallow  water  indicates 
all  water  depths  for  which  sound  Interacts  significantly  with  the 
ocean  bottom.  The  separation  frequency  of  500  Hz  between  the  low- 
and  high-frequency  regimes  is  arbitrarily  chosen. 

When  indicating  the  applicability  of  a  propagation  model  to  a 
given  type  of  environment  we  take  into  consideration  limitations 
in  the  underlying  theory.  Ray  models  are  applicable  mainly  to 
«.  high-frequency  propagation.  Only  ray  and  PE  models  accurately 
handle  a  range-dependent  environment.  The  normal-mode  model 
treats  range  dependence  in  the  adiabatic  approximation.  When 
indicating  a  model's  practicality  we  consider  exclusively  the  com¬ 
putation  time,  which,  of  course,  depends  on  the  required  accuracy. 
The  computation  time  increases  with  both  frequency  and  water  depth 
for  wave  models  (mode,  FFP,  PE),  while  the  time  is  relatively 
independent  of  these  parameters  for  the  ray  models.  Likewise, 
computation  time  is  proportional  to  the  number  of  profiles  in  a 
range-dependent  environment  for  both  ray  and  mode  models,  while  a 
PE  model  takes  essentially  the  same  time  for  range-dependent  and 
range-independent  environments. 

Full-box  shading  in  Fig,  9  means  that  a  model  is  applicable  as 
well  as  practical.  On  the  other  hand,  if  a  box  is  only  partially 
shaded,  it  means  that  the  model  is  applicable  with  caution 
(theoretical  limitations),  or  that  computation  times  are 
excessive.  The  above  judgements  are,  of  course,  relative.  For 
instance,  in  our  first  evaluation  some  columns  contained  no  appli¬ 
cable  models.  In  these  columns  we  therefore  selected  the  model  we 
felt  was  the  most  practical  end  denoted  it  by  a  fully  shaded  box. 
For  a  column  where  more  than  one  ‘  box  is  fully  shaded,  the  choice 
of  model  will  depend  on  the  actual  models  on  hand,  the  running 
time,  input/output  options  available,  etc. 

Since  the  various  models  are  approximate  solutions  of  the  wave 
equation,  it  is  valuable  to  check  the  validity  of  these  approxima¬ 
tions  by  doing  an  inter-model  comparison  for  situations  where  all 
four  models  are  considered  applicable.  Returning  to  Fig.  9,  we 
note  that  all  models  should  handle  ?  range-independent  shallow- 
water  environment  at  around  500  Hz,  even  though  the  mode  and  FFP 
models  are  designated  most  applicable  (fastest). 

An  example  of  a  transmission  loss  computation  for  shallow 
water  is  given  in  Fig.  10.  Here  an  isovelocity  water  column 
(1500  m/s)  is  100  m  deep  and  both  source  (SD)  and  receiver  (RD) 
are  at  mid-depth.  The  bottom  is  considered  homogeneous,  with  a 
sound  speed  of  1550  m/s,  a  density  of  1.2  g/cm3,  and  an  atte¬ 
nuation  of  1  dB/wavelength.  We  see  from  Fig.  10  that  the  three 
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Fig.  lO  Inter-model  comparison  for  shallow-water  environment 


wave  models  (mode,  FFP,  PE)  give  virtualfy  identical  results  both 
for  level  and  for  the  multipath  interference  structure  as  a  func¬ 
tion  of  range.  The  ray  model ,  though  it  cannot  reproduce  the 
interference  pattern,  does  yield  the  same  approximate  level 
(dashed  line). 

Many  consistency  tests  of  acoustic  models  have  to  be  carried 
out  in  order  to  check  the  complex  computer'  programs  within  the 
framework  of  the  theoretical  limitations  particular  to  each  model 
<39>.  A  positive  outcome  of  an  inter-model  comparison  helps  us  to 
gain  confidence  in  particular  numerical  models.  However,  we 
should  remember  that  the  final  check  on  an  acoustic  model  is  a 
comparison  with  experimental  data.  This  demonstrates  whether  or 
not  the  model  includes  all  the  physics  necessary  for  explaining 
and  understanding  sound  propagation  in  a  real  ocean  environment. 
A  series  of  model/data  comparisons  can  be  found  in  <40-43>. 
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5  SPECIAL  MODEL  APPLICATIONS 

This  section  is  dedicated  to  a  study  of  some  general  wave- 
propagation  problems,  for  which  Illustrative  numerical  solutions 
can  be  obtained  in  a  straightforward  manner  using  one  of  the 
aforementioned  ocean-acoustic  models.  We  shall  first  address  the 
problem  of  reflection  of  a  gaussian  beam  of  arbitrary  width  at  a 
fluid/solid  interface  near  the  Rayleigh  angle,  where  a  leaky  sur¬ 
face  wave  is  excited  causing  a  complex  reflection  pattern  with 
beam  splitting  and  beam  displacement.  This  is  a  well-researched 
problem  in  both  optics  and  acoustics,  which  we  can  easily  solve 
with  the  fast-field  program.  Next  we  study  propagation  in  range- 
dependent  waveguides  using  the  parabolic  equation  technique.  We 
shall  address  acoustic  radiation  from  a  duct  into  free  space  as 
well  as  mode  coupling  in  tapered  waveguides. 

5.1  Beam  reflection  at  fluid/solid  interface 


Bounded  beam  reflection  near  a  critical  angle  is  a  subject 
that  has  received  considerable  attention  in  the  past,  both  within 
the  field  of  electromagnetics  <44-48>  and  acoustics  <49-55>.  The 
phenomenon  of  interest  has  mainly  been  the  displaced  reflected 
beam,  while  the  transmitted  field  has  been  studied  in  much  less 
detail.  There  are  basically  two  different  wave  phenomena  that  can 
account  for  the  observed  features  of  the  reflected  field.  One  is 
the  excitation  of  a  lateral  wave  at  the  interface  when  a  beam  of 
finite  width  is  incident  on  the  interface  at  grazing  angles  lower 
than  the  critical  angle.  The  reflected  field  is  then  composed  of 
contributions  from  both  the  specular  reflection  and  the  lateral 
wave  field,  causing  an  apparent  lateral  displacement  of  the 
reflected  beam.  In  acoustics,  the  lateral  wave  phenomenon  is 
associated  with  beam  reflection  at  fluid/fluid  interfaces. 

The  second  phenomenon  is  associated  with  the  excitation  of  a 
leaky  surface  wave,  which  again  complicates  the  reflection  pattern 
when  added  to  the  specularly  reflected  field.  This  phenomenon 
occurs  in  acoustics  when  a  bounded  beam  is  incident  on  a 
fluid/solid  interface  just  below  the  shear  critical  angle. 

Various  aspects  of  bounded  beam  physics  have  been  studied 
theoretically  as  well  as  experimentally  in  the  past  35  years.  The 
lateral  displacement  of  a  reflected  light  beam  was  first  observed 
by  Goos  and  Hanchen  <44>,  and  the  phenomenon  is  therefore  often 
referred  to  as  the  Goos-Hanchen  effect.  Several  theoretical 
papers  have  addressed  the  beam  reflection  problem,  for  instance 
<45-49>,  though  always  with  some  theoretical  limitations,  such  as 
bearawidth  large  compared  to  the  wavelength,  lossless  media, 
parallel  beams  of  particular  shape  (gaussian),  etc. 

We  shall  here  apply  the  fast-fi^ld  program  (FFP)  to  numeri¬ 
cally  solve  the  reflection  problem,  and,  as  we  shall  see,  this 
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technique  treats  the  complete  problem  without  any  of  the  above 
limitations.  We  are  going  to  study  reflection  at  a  fluid/solid 
interface,  where  the  reflection  phenomenon  is  associated  with  the 
excitation  of  a  leaky  surface  wave  (Rayleigh  wave).  However,  the 
FFP  model  could  as  well  be  applied  to  reflection  and  transmission 
for  a  solid  plate  In  a  liquid  <52— 53>  or  to  the  transmission 
through  a  fluid/fluid  boundary  <54-55>. 

The  FFP  model  <13>  provides  an  exact  numerical  solution  for 
the  acoustic  field  generated  by  a  point  source  in  a  multilayered 
fluid/solid  environment  (Sect.  3.1).  Wave  attenuation  for  both 
compressional  and  shear  waves  is  included  in  the  theory.  We  have 
modified  the  standard  code  to  efficiently  solve  the  system  of 
equations  for  a  number  of  point  sources,  equldistantly  spaced  and 
forming  a  vertical  line  array.  The  total  acoustic  field  is  found 
by  superposition  of  the  contributions  from  individual  point  sour¬ 
ces.  Hence,  in  this  model  the  beam  is  generated  by  a  vertical 
soiree  array,  and  the  beam  direction  is  varied  by  appropriately 
phasing  the  source  elements.  By  using  a  gaussian  amplitude 
weighting  across  the  array,  a  gaussian  beam  can  be  generated,  pro¬ 
pagating  at  any  angle  with  respect  to  the  horizontal.  By  varying 
the  array  distance  from  the  interface  and  the  number  of  source 
elements  (half-wavelength  spacing) ,  a  beam  of  arbitrary  width  can 
be  generated,  and  we  can  obtain  parallel,  diverging,  or  converging 
beams,  as  we  wish.  Hence,  the  model  is  very  general  in  concept, 
and  should  handle  the  reflection  problem  accurately  for  any  beam- 
width  and  for  multilayered  systems  (plates)  as  well  as  the  reflec¬ 
tion  and  transmission  at  a  single  fluid/solid  interface.  The 
solution  is  for  a  plane  geometry,  l.e.  a  two-dimensional  beam. 

We  first  consider  computational  results  for  a  water/steel 
interface.  Information  on  material  parameters  were  taken  from 
Breazeale  et  al  <51>,  and  the  FFP  results  will  be  compared  with 
the  theory  of  Bertoni  and  Tamlr  <47>,  here  named  the  BT  theory. 
Results  for  reflection  at  the  Rayleigh  angle  (59.35°)  are  given  in 
Fig.  11  for  two  different  beamwidths. 

Computations  were  done  at  20  kHz,  and  the  acoustic  field  is 
displayed  as  iso-loss  contours  given  in  arbitrary  decibels  (low 
values  correspond  to  high  intensity).  We  are  considering  incoming 
parallel  gaussian  beams,  where  the  beamwidth  is  10  and  64  wave¬ 
lengths  respectively.  Note  in  Fig.  11a  that  the  narrow  beam,  when 
reflected,  is  being  split  up  in  two  beams  separated  by  a  "null 
strip”.  The  beam  to  the  left  is  displaced  slightly  backwards, 
while  the  rightmost  beam  is  displaced  forwards  18  wavelengths. 
This  behavior  is  in  complete  agreement  with  the  BT  theory.  With 
increasing  beamwidth,  the  energy  shifts  to  the  forwardly  displaced 
beam,  and  we  end  up  (Fig.  lib)  with  a  single-beam  reflection  with 
a  displacement  of  50  wavelengths.  The  general  behavior  seen  in 
these  examples  is  predicted  by  Bertoni  and  Tamlr,  and  the  measured 
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Fig.  11  Beam  reflection  at  water/steel  interface  at  the 
Rayleigh  angle. 

a)  narrow  beam  of  width  lo\ , 

b)  broad  beam  of  width  64\. 

\ 

displacement  versus  beamwidth  gives  data  points  that  fall  exactly 
on  the  displacement  curve  computed  from  the  BT  theory. 


A  more  challenging  Investigation  was  to  explain  the 
disagreement  between  theory  (BT)  and  experimental  results  for  alu¬ 
minium  oxide,  as  reported  by  Breazeale  et  al  <51>.  We  proceed  to 
calculate  the  Raylelgh-wave  properties ' for  a  water/alumlnlum-oxlde 
interface  using  the  material  parameters  given  in  <51>.  For  a 
water  speed  of  1490  m/s,  a  frequency  of  2  MHz,  and  with  realistic 
attenuation  coefficients  in  both  media,  we  found  the  leaky  surface 
wave  to  have  a  phase  velocity  of  5825. A  m/s  corresponding  to  an 
angle  of  75.18®.  The  leakiness  of  the  surface  wave  is  given  by 
the  imaginary  part  of  the  wave  number,  which  Is  calculated  to  be 
23.3  a".  These  Rayleigh-wave  properties  were  determined  from  a 
separate  computer  code  that  finds  the  poles  for  the  complex 
reflection  coefficient. 
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The  next  step  was  to  generate  a  parallel  beam  (plane  wave 
front)  with  a  half-width  of  16.6  tin  as  used  in  the  experiment. 
The  computed  field  at  the  Rayleigh  angle  Is  shown  in  Fig.  12. 
Again  we  see  that  the  reflected  beam  is  spilt  In  two,  and  we  also 
notice  the  leading  radiated  field  associated  with  the  surface 
wave.  In  fact,  full  Information  about  the  Rayleigh  wave  can  be 
read  off  directly  from  the  rightmost  contour  lines  (46  to  54  dB). 
They  have  a  slope  of  75.2°  with  the  horizontal  (the  Rayleigh 
angle),  and  the  field  decay  is  23.3  Nepers/m  (the  imaginary  part 
of  the  wavenumber).  We  found  that  this  property  of  the  reflected 
field  is  independent  of  the  angle  of  incidence  of  the  incoming 

befm,  and,  as  we  shall  see,  it  is  a  valid  criterion  for  deter¬ 
mining  the  Rayleigh-wave  properties  also  for  a  diverging  beam. 

The  computed  intensity  distribution  in  the  reflected  beam  is 
given  in  Fig.  13,  normalized  with  the  amplitude  of  the  incoming 

beam  measured  at  the  interface  (dotted  profile  at  75.2°).  The 
vertical  dashed  line  corresponds  to  the  position  of  a  specular 
reflection.  The  results  are  computed  40  cm  above  the  interface, 
corresponding  to  a  horizontal  cut  through  the  contour  plot 
(outside  the  frame).  Hence,  this  display  differs  slightly  from 
the  experimental  results,  which  were  obtained  measuring  perpen¬ 
dicular  to  the  beam  direction.  This  difference  is  however  minor 
at  these  angles.  We  see  from  Fig.  13  the  expected  behavior,  i.e. 
a  single  specularly  reflected  beam  when  moving  1.5  to  2.0°  away 

from  the  Rayleigh  angle.  We  also  notice  that  the  interference 

null  is  strongest  at  the  Rayleigh  angle  (75.2°)  as  predicted  by 
the  BT  theory  for  parallel  beams.  Hence  these  results  confirm 
that  the  Rayleigh  angle  for  the  chosen  material  parameters  is 
75.2°. 


Fig .  12 

Reflection  of  parallel 
beam  at  water/aluminium- 
oxide  interface 
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We  now  consider  a  slightly  diverging  beam  (curved  wave  front) 
as  actually  used  in  the  experiment.  The  beamwidth  is  again 
16.6  ram,  and  the  divergence  measured  at  the  3  dB  down-points  is 
2.7°  at  the  Interface.  Our  results  for  the  reflected  field  at  the 
Rayleigh  angle  is  given  in  Fig.  14.  We  notice  that  the  reflected 
field  is  now  much  more  complicated,  with  essentially  3  reflection 
peaks.  However,  the  undisturbed  leading  edge  of  the  radiated  sur¬ 
face  wave  exhibits  the  same  properties  as  before,  i.e.  the  iso- 
intensity  contours  have  a  slope  of  75.2°  with  the  horizontal,  and 
the  field  decay  is  23.3  Nepers/m  parallel  to  the  interface. 

The  computed  field  intensity  40  cm  above  the  interface  is 
given  in  Fig.  15  as  a  function  of  the  incoming  beam  angle.  We 
notice  the  many  interference  lobes  now  present  in  the  reflected 
beam.  In  this  case  the  reflection  pattern  at  the  Rayleigh  angle 
(75.2°)  has  no  particular  features,  such  as  a  pronounced  inter¬ 
ference  null,  that  makes  it  possible  to  easily  determine  the 
Rayleigh  angle.  In  fact,  when  searching  for  the  strongest  inter¬ 
ference  null,  one  finds  this  to  occur  at  approximately  75.8°, 
which  is  very  close  to  the  angle  (7  5.7°)  designated  the  Rayleigh 
angle  in  the  experiment  <5 1 > ,  using  the  above  (wrong)  criterion. 
Hence,  we  may  conclude  from  this  set  of  curves,  that  an  experimen¬ 
tal  verification  of  the  Rayleigh  angle  can  most  easily  be  done  by 
measuring  the  slope  of  the  lsoloss  contours  in  the  undisturbed 
leaky  wave  field  as  seen  an  the  contour  plots. 

A  detailed  comparison  between  our  theoretical  result  at  75.8° 
and  the  experimental  result  given  in  Ref.  <51>,  shows  a  clear 
improvement  over  the  prediction  obtained  from  the  BT  theory. 
However,  some  disagreement  on  peak  levels  still  exists,  which  we 


Fig.  14 

Reflection  of  diverging 
beam  at  water/aluminium- 
oxide  interface 
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Fig.  15  Reflectivity  pattern  versus  grazing  angle  for  diverging 
beam  incident  on  a  water/aluminium-oxide  interface 
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speculate  could  be  due  to  Inaccurate  knowledge  of  the  material 
parameters  for  aluminium  oxide. 

The  beam  reflection  calculations  with  the  FFP  model  is  an 
ongoing  project.  We  are  mainly  interested  in  studying  the  reflec¬ 
tion  process  for  very  narrow  (focused)  beams,  where  the  reflection 
process  is  strongly  influenced  by  diffraction  effects. 

5.2  Propagation  from  duct  into  free  space 

The  study  of  sound  propagation  in  a  range-dependent  environ¬ 
ment  is  a  fascinating  subject,  to  which  we  will  devote  the 
remainder  of  this  paper.  We  have  chosen  to  solve  a  series  of 
relatively  simple  propagation  situations  suitable  for  the  parabo¬ 
lic  equation  technique.  As  pointed  out  in  Sect.  3.3,  the  PE 
method  is  limited  to  propagation  within  ±  20°  with  respect  to  the 
horizontal,  and  back-scattering  is  neglected  in  the  solution.  The 
PE  calculations  have  been  done  for  a  cylindrical  geometry. 

One  of  the  simplest  range-dependent  problems  in  acoustics  is 
the  radiation  of  sound  from  a  symmetric  duct  into  free  space.  PE 
solutions  to  this  problem  are  shown  in  Fig.  16.  The  duct  is  100  m 
wide  with  a  sound  speed  of  1500  m/s  and  a  density  of  1  g/cm^.  The 
infinitely  thick  duct  wall  has  a  speed  of  1550  m/s.  There  is  no 
density  change  in  the  problem,  and  material  losses  are  neglected. 
At  a  frequency  of  50  Hz,  there  are  two  propagation  modes  in  the 
duct. 


Figure  16a  shows  the  computed  field  from  a  source  placed  in 
the  middle  of  the  duct.  In  this  case  only  the  symmetric  1st  mode 
is  excited,  radiating  symmetrically  into  free  space  beyond  a  range 
of  2  km.  The  initial  modal  field  for,  the  PE  calculation  was 
supplied  by  a  normal-mode  model.  Next  we  moved  the  source  to  a 
depth  of  27  m  below  the  center  of  the  'duct,  which  gives  equal 
excitation  of  the  two  modes.  Figure  16b  shows  that  the  radiated 
field  is  now  split  up  into  two  almost  symmetric  beams.  We  would 
expect  the  radiated  field  to  be  determined  by  the  field  distribu¬ 
tion  across  the  duct  opening.  This  is  confirmed  in  Fig.  16c, 
where  the  duct  has  been  truncated  at  1.5  km  range;  we  now  obtain 
an  asymmetric  radiation  pattern  with  most  of  the  energy  being  con¬ 
tained  in  the  down-going  beam. 

The  above  numerical  results  provide  considerable  physical 
Insight  into  the  duct  radiation  problem.  The  examples  were  chosen 
so  as  to  facilitate  a  physical  interpretation  of  the  contour 
plots.  However,  more  complex  situations  could  easily  be  Investi¬ 
gated. 
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5.3  Up-slope  propagation 

We  now  proceed  to  study  propagation  over  sloping  bottoms  at  a 
sufficiently  low  frequency  that  phenomena  such  as  mode  cutoff  and 
mode  conversion  can  be  Investigated  in  some  detail.  First  we  con¬ 
sider  up-slope  propagation  for  the  environment  given  in  Fig.  17. 
The  water/bottom  interface  is  indicated  on  the  contour  plot  by  the 
heavy  line  starting  at  350  m  depth  and  moving  towards  the  surface 
beyond  a  range  of  10  km.  The  bottom  slope  is  0.85°.  The  fre¬ 
quency  is  25  Hz  and  the  source  depth  is  150  m.  The  water  is  taken 
to  be  Isovelocity  with  a  speed  of  1500  m/s,  while  the  bottom  is 
characterized  by  a  speed  of  1600  m/s,  a  density  of  1.5  g/cm^,  and 
an-oattenuation  of  0.2  dB/wavelength. 

Before  interpreting  the  contour  plot,  let  us  have  a  look  at 
the  simplified  sketch  in  the  upper  part  of  Fig.  17.  Using  the 
ray/mode  analogy,  a  given  mode  can  be  associated  with  up-  and 
down-going  rays  with  a  specific  grazing  angle  corresponding  to  a 
given  mode.  As  sound  propagates  up  the  slope,  the  grazing  angle 
for  that  particular  ray  (mode)  increases,  and  at  a  certain  point 
in  range  the  angle  exceeds  the  critical  angle  at  the  bottom, 
meaning  that  the  reflection  loss  becomes  very  large  and  that  the 
ray  essentially  leaves  the  water  column  and  starts  propagating  in 
the  bottom.  The  point  in  range  where  this  happens  corresponds  to 
the  cutoff  depth  for  the  equivalent  mode. 
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Fig.  17  Up-slope  propagation  showing  discrete  mode  cutoffs 
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To  emphasize  the  main  features  in  Fig.  17 ,  we  have  chosen  to 
display  contour  levels  between  70  and  100  dB  in  2  dB  intervals. 
Thus  high-intensity  regions  (loss  <  70  dB)  are  given  as  blank 
areas  within  the  wedge,  while  low-intensity  regions  (loss 
>  100  dB)  are  given  as  blank  areas  in  the  bottom.  The  PE  solution 
was  here  started  off  by  a  gaussian  initial  field,  and  there  are 
four  propagating  modes.  The  high  intensity  in  the  bottom  at  short 
ranges  (<  10  km)  corresponds  to  the  radiation  of  "continuous 
modes”  into  the  bottom.  As  sound  propagates  up  the  slope  we  see 
four  well-defined  beams  in  the  bottom,  one  corresponding  to  each 
of  the  four  modes.  This  phenomenon  of  energy  leaking  out  of  the 
■*-  propagation  channel  as  discrete  beams  has  been  confirmed  experi¬ 
mentally  <56>,  and  a  detailed  study  of  this  phenomenon  using  the 
PE  method  has  been  reported  elsewhere  <57>. 

This  particular  example  of  mode  coupling  where  discrete  modes 
trapped  in  the  water  column  couple  into  "continuous"  modes  propa¬ 
gating  in  the  bottom,  has  received  much  attenation  recently; 
several  theoretical  papers  <3l>  and  <58-60>  have  appeared  offering 
solutions  to  the  wedge-propagation  problem. 

We  now  proceed  to  study  the  problem  of  mode  coupling  within 
the  wedge  itself.  Strong  coupling  can  be  achieved  either  by 
Increasing  the  bottom  slope  or  by  increasing  the  frequency.  It  is 
the  latter  case  that  will  be  considered  here.  Figure  18  shows  PE 
calculations  for  a  simple  wedge  problem.  The  initial  water  depth 
is  100  m  and  the  bottom  slope  is  2°.  The  water  is  taken  to  be 
isovelocity  with  a  speed  of  1500  m/s,  while  the  bottom  has  a  speed 
of  1550  m/s.  There  is  no  density  change  in  the  problem,  and 
material  losses  are  neglected.  The  source  depth  is  50  m,  and  the 
initial  field  for  the  PE  calculation  was  supplied  by  a  normal-mode 
model.  The  two  plots  are  for  source  frequencies  of  50  and  500  Hz, 
respectively.  In  both  cases,  only  the  first  mode  was  propagated 
up  the  slope.  , 

We  notice  in  Fig.  18a  that  no  mode  coupling  takes  place  at  a 
frequency  of  50  Hz.  The  contour  lines  are  smooth,  indicating  that 
the  local  mode  at  range  zero  adapts  to  the  changing  water  depth 
until  it  reaches  the  cutoff  depth,  where  the  energy  radiates  into 
the  bottom.  This  is  an  example  where  propagation  is  well 
described  by  adiabatic  mode  theory  <2 1 — 26> . 

Strong  mode  coupling  occurs  when  we  increase  the  frequency 
to  500  Hz  (Fig.  18b).  The  field  within  the  wedge  becomes  compli¬ 
cated,  and  so  does  the  radiation  pattern  into  the  bottom.  We  have 
done  no  attempt  to  analyse  this  complex  mode-coupling  problem  in 
detail,  but  the  PE  technique  could  certainly  be  used  for  such  a 
study. 
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Fig.  18  Up-slope  propagation  over  a  constant  2°  slope 

a)  50  Hz,  no  mode  coupling, 

b)  500  Hz,  strong  mode  coupling . 


5.4  Down-slope  propagation 

We  now  consider  the  problem  of  down-slope  propagation  as 
illustrated  in  Fig.  19.  The  initial  water  depth  is  50  m  and  the 
bottom  slope  is  5°.  The  water  column  is  isovelocity  with  a  speed 
of  1500  m/s,  while  the  bottom  has  a  speed  of  1600  m/s.  The  den¬ 
sity  ratio  between  bottom  and  water  is  1.5,  and  a  wave  attenuation 
of  0.5  dB/wavelength  has  been  included  in  the  bottom.  The  two 
contour  plots  are  for  source  frequencies  of  25  and  500  Hz,  respec¬ 
tively.  The  Initial  fields  for  the  PE  calculations  were  supplied 
by  a  normal-mode  model,  and,  in  both  cases,  <nly  the  first  mode 
was  propagated  down  the  slope. 
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An  analytical  solution  of  Eq.  19  is  of  course  possible 
leading  to  closed-form  expressions  for  the  arbitrary  functions; 
but  for  more  than  a  few  layers  this  procedure  would  be  incon¬ 
venient.  Further,  the  Hankel  transforms  do  not  lead  to  closed- 
form  solutions,  but  need  numerical  evaluation.  Thus  the  most 
general  way  to  proceed  is  to  create  a  numerical  model  based 
directly  on  the  system  of  equations  (Eq.  19).  Such  a  model  is 
described  in  the  next  chapter. 

2  THE  NUMERICAL  MDDEL  SAFARI 

The  numerical  evaluation  of  the  Hankel  transform  necessitates 
a  truncation  and  a  discretization  in  the  horizontal  wavenumber  s. 
As  can  be  observed  from  Eqs.  8  and  17  ,  the  source  terms  decay  expo'- 
nentially  for  s  going  towards  infinity.  As  the  source  terms  form 
the  righthand  side  of  Eq.  19  the  arbitrary  functions  will  behave  in 
the  same  way.  It  is  therefore  possible  to  truncate  the  integration 
interval  in  accordance  with  any  accuracy  demands.  The  fast-field 
technique  introduced  by  Marsh  <7>  can  then  be  used  to  evaluate  the 
Hankel  transforms. 

The  Bessel  functions  are  expressed  in  terms  of  Hankel  func¬ 
tions 


Jm(rs)  -  1/2  (Hd)(rs)  +  H(2)(rs)) 
m  m 


(20) 


and  each  integral  is  split  into  two.  As  only  outgoing  waves  are 
considered,  the  integrals  involving  .  H^l(rs)  are  neglected,  and 

H^2^(r8>  is  replaced  by  its  asymptotic  form 


H(2)(r8)  ~ 

m  rs+— 


(21) 


The  integration  over  the  truncated  interval  can  now  be  per¬ 
formed  by  means  of  the  fast  Fourier  transform,  and  the  actual  field 
parameter  is  found  at  a  number  of  ranges  equal  to  the  number  of 
discrete  wavenumbers  considered. 

DiNapoli  and  Deavenport  <3>  have  compared  Marsh's  method  to 
the  technique  introduced  by  Tsang,  Brown,  Kiang  and  Simmons  <8>, 
which  does  not  use  the  asymptotic  form  of  the  Hankel  function. 
Significant  differences  were  found  only  for  very  short  ranges. 

The  kernels  in  the  Hankel  transforms  are  now  needed  only  for 
a  limited  number  of  discrete  values  of  s.  Kutschale  <4>  used  a 
Green 's-function  approach  based  on  the  Thomson-Haskell  matrix 


( 


49 

3SN3dX3  1N3INNH  3AOD  IV  M  vm<ion,»  >>| 


SACLANTCEN  SR-83 


arz<r.z>  |  n  =  Un  /  {2san(A“ 


_zan  A+  zar 
e  -  A  e 

n  n 


-(2s^-k^)(B_  e  n  +  e  n)}  s  J.(rs)  ds. 


n  n 


(15) 


In  the  case  of  a  fluid  layer  the  shear  stiffness  pn  vanishes, 
an4  only  the  potential  is  involved.  The  displacements  follow 
directly  from  Eqs.  12  and  13  by  setting  B"  and  B+  to  zero.  The 
shear  stress  is  identically  zero,  whereas  Eq .  14  has  to  be 
replaced  by 


(r»z)  n  =  -^n^  /  (a  e  +  A+  e  n}  s  T0(rs)ds. 


’zz 


(16) 


The  source  is  assumed  to  be  in  layer  number  m  at  depth  zs. 
In  the  absence  of  boundaries,  the  field  produced  in  layer  m  would 
be  <7>: 


iSu  .  "  I  2  z8  I  am 

#s(r,z)  =  — —  /  - - -  s  JQ(rs)ds  , 


4* 


Ag(r,z)  =  0  , 


(17) 

(18) 


where  is  the  source  strength.  If  Eq.  17  Is  inserted  into 

Eq.  1,  expressions  similar  to  Eqs.  12  and  13  are  obtained  for  the 
displacements,  and  again  Hooke’s  law  yields  expressions  like  Eqs. 
14  and  15  for  the  stresses  involved  in  the  boundary  conditions. 

For  each  value  of  the  range,  r,  the  boundary  conditions  must 
be  satisfied.  In  the  upper  and  lower  half-spaces  the  arbitrary 
functions,  with  superscript  -  and  +  respectively,  must  vanish  due 
to  the  radiation  condition.  At  each  interface,  w  and  ozz  must  be 
continuous,  and  at  all  solid/liquid  interfaces  the  shear  stresses 
must  vanish.  At  solid/solid  interfaces,  w,  u,  ozz  and  arz  must 
be  continuous.  This  yields  a  linear  system  of  equations  in  the 
arbitrary  functions,  to  be  satisfied  at  each  horizontal 
wavenumber: 


Cij(8)  •  Aj(s)  *  Ri(s).  (19) 

The  vector  Aj(s)  contains  all  the  non-vanishing  arbitrary 
functions,  Cij(s)  is  the  coefficient  matrix,  and  Rj(s)  contains 
the  contributions  from  the  source.  When  the  arbitrary  functions 
are  found,  the  field  parameters  at  any  depth  and  ranee  can  be 
obtained  from  the  Hankel  transforms  (Eqs.  12  to  16)  plus  the 
source  contributions  (if  the  source  and  receiver  ere  t-  th<-  same 
layer). 
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an(s),  Bn(s)  are  defined  as 


*n(s) 


*«<■> 


/ (s2-h2) , 

s  >  Re( hn| 

i  /(h2-s2), 

s  <  Re(hn} 

/ (s2-kn) , 

s  >  Re{kn} 

i  /(k^-s2). 

s  <  Re{kn} 

(8) 


(9) 


The  wavenumbers  hn  and  kn  for  the  compressional  and  shear 
waves,  respectively,  are  defined  by 


“2Pn 

^n+2M! 

M 

CM 

(-s-)2 

u2pn 

m  — —  — 

n 

Orn 

Hn 

If 

Eqs.  6 

and  7 

(10) 


(11) 


expressions  are  obtained  for  the  particle  displacements: 


w(r,e)  |  -  /  {-on  A”  e  Z°n  +  on  & 

n  n  n 


zar 


+  8B“  e  +  s  B+  e  ^n}  s  J0(rs)ds,  (12) 


-  -“Wn  -  s  A+  •“» 


u(r,z)  |n"/  {-s  A-  e 
o  n 

+  Bn  B"  e  n  -  8n  B+  e  n|  s  Ji(rs)ds.  (13) 


The  stress  components  involved  in  the  boundary  conditions 
follow  from  Hooke's  law: 


'zz 


(r,z)  |  n  -  un  /  {(2s2-k2)(A“  e  Z  n  +  A+  e  n) 

Z^n  +  B+  e  n){  s  J0(rn)ds, 


+  2sBn(-B"  e 
n 


(14) 
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*n  ~  0 
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(?2 
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o. 
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in  which  Cl  and  Gj.  are  the 
shear  waves,  respectively: 

velocities  of 

the  compress! onal  and 

r>2 

^n  +  2un 

(4) 

CLn 

Pn 

CT 

Tn 

Pn 

* 

(5) 

In  the  present  case  the  field  is  asymmetric  due  to  the  posi¬ 
tioning  of  the  source  on  the  axis,  and  the  angular  displacement  v 
vanishes  everywhere.  It  is  then  clear  from  Eq.  1  that  the  poten¬ 
tial  7n  must  be  constant  and  can  be  excluded. 

In  the  following,  only  vibrations  with  angular  frequency  w 
will  be  considered;  displacements,  stresses  and  potentials  can 
then  be  expressed  In  complex  form  with  the  common  factor  exp(iwt). 
This  factor  will  not  be  included  in  '  the  following.  The 
viscoelastic  damping  can  now  be  accounted  *for  by  allowing  the  Lame 
constant  to  be  complex.  After  use  of  the  Hankel  transform  on  the 
wave  equations,  the  following  integral  representations  are 
obtained  for  the  potentials: 

*n(r,z)  -  /  {A“(s)e  Z<ln^  ^  +  A+(s)e  JQ(rs)s  ds  (6) 


An(r,z)  *  /  {B“(s)e  ^  +  B+(s)e  ^  JQ(rs)s  ds,  (7) 

o  n  n 

where 

Jm  is  the  Bessel  function  of  the  first  kind  and  order  m, 

A”,  A+,  B~  and  B+  are  arbitrary  functions  in  the  horizontal  wave- 
n  n  n  n 

number  S. 
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Even  with  one  source/recetver  combination,  the  calculation 
speed  has  been  Improved  by  up  to  one  order  of  magnitude.  This 
gain  In  speed  and  performance  has  made  the  FFP-technique  a  more 
realistic  alternative  to  other  numerical  techniques,  and  in  areas 
where  this  technique  is  the  only  possibility,  the  range  of 
solvable  problems  has  been  increased  considerably. 

In  the  following  the  model  and  its  mathematical  background  is 
described,  and  two  examples  of  its  use  are  gtven.  First  it  is 
used  to  clarify  some  of  the  basic  properties  of  the  slow  seismic 
interface  waves  that  can  be  observed  in  shallow  water  environ¬ 
ments.  Synthetic  seismograms  are  produced  and  compared  qualitati¬ 
vely  to  experimental  data.  Then  the  model  is  used  to  determine 
the  reflection  of  a  pulsed  beam  from  a  water/aluminium-oxide 
interface. 


1  THE  MATHEMATICAL  MODEL 


The  mathematical  model  is  based  on  the  assumption  that  the 
water  column  and  the  bottom  consist  of  a  series  of  range- 
independent  layers.  All  materials  are  considered  to  be  homoge¬ 
neous  and  isotropic  elastic  continue  with  Lame  constants  An  and  pn 
and  density  pn.  The  subscript  refers  to  layer  number  n.  The 
damping  mechanisms  are  assumed  to  be  linear  viscoelastic. 


A  cylindrical  coordinate  system  {r,0,z}  is  introduced  with 
the  z-axis  going  through  the  source  and  being  positive  downwards 
(Fig.  1).  The  representation  of  the'  cylindrical  displacement  com¬ 
ponents  (u,v,w)  in  terms  of  scalar  potentials  and  the  subsequent 
expression  of  these  as  Hankel  transforms  closely  follows  the  pre¬ 
sentation  given  by  Schmidt  and  Krenk  <6>;  hence,  only  an  outline 
will  be  given  here.  If  body  forces  are  neglected,  the  displace¬ 
ment  equation  of  motion  will  be  satisfied  if  the  displacement  com¬ 
ponents  in  layer  n  are  expressed  in  terms  of  three  scalar 
potentials  {*n»*n»^n}  88 
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where  the  potentials  satisfy  the  wave  equations: 
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on  the  geometry  (plane  or  cylindrical)  the  field  Is  decomposed 
into  plane  or  conical  waves.  The  field  corresponding  to  each 
horizontal  wavenumber  Is  then  found  by  matching  of  boundary  con¬ 
ditions  at  each  Interface,  and  the  total  field  is  found  by  super¬ 
position.  In  the  case  of  cylindrical  geometry  the  superposition 
is  given  by  a  Hankel  transform  integral.  By  replacing  the  Hankel 
function  with  its  large  argument  approximation,  this  is  changed  to 
a  ''Fourier  integral,  which  after  truncation  can  be  evaluated  by 
means  of  a  fast  Fourier  technique.  This  approximation  of  the 
Hankel  transform  was  originally  introduced  by  Marsh  <2>  and  is 
usually  called  the  fast  field  technique. 

The  major  part  of  the  numerical  effort  is  related  to  the 
solution  of  the  depth-separated  wave  equation.  DiNapoli  <3> 
introduced  an  FFP-code  that  performed  the  solution  very  effi¬ 
ciently  by  means  of  recurrence  relations  for  the  hypergeometric 
functions.  However,  his  approach  allows  only  for  fluid  layers, 
and  in  such  cases  other  models  will  often  be  more  convenient,  e.g. 
normal-mode  models. 

The  main  advantage  of  FFP-models  are  their  ability  to  treat 
problems  involving  solid  layers  <1>.  The  first  model  to  include 
this  feature  was  introduced  by  Kutschale  <4>.  As  shear  properties 
are  very  important  for  sound  propagation  in  shallow  water  <5> , 
this  model  has  been  used  extensively  for  such  problems  during  the 
last  decade.  Kutschale  solved  the  depth-separated  wave  equation 
using  a  Green ' s-function  approach  based  on  the  Thomson-Haskell 
matrix  method.  However,  this  method  allows  for  only  one  source 
and  receiver  at  a  time,  and  even  in  such  cases  the  computations 
are  rather  extensive,  often  the  modelling  of  pulse  propagation 
becomes  impractical  since  this  has  to  be  carried  out  by  means  of 
Fourier  synthesis  involving  many  frequency  components. 

Here  a  new  FFP-code,  called  SAFARI,  is  introduced.  Instead 
of  using  the  Thomson-Haskell  method,  a  system  of  linear  equations 
in  the  unknown  potentials  is  set  up  and  solved  at  each  horizontal 
wavenumber.  When  the  layer  series  is  given;  it  is  possible  from 
the  boundary  conditions  to  determine  a  priori  a  mapping  between 
the  equations  to  be  satisfied  at  each  Interface  and  a  global  set 
of  equations.  A  significant  number  of  wavenumber-independent 
expressions  can  then  be  evaluated  once.  The  mapping  technique  is 
very  similar  to  the  one  used  in  finite-element  programs. 

The  source  contributions  are  mapped  into  the  righthand  side 
of  the  global  equations,  and  there  is  no  theoretical  limit  for  the 
number  of  sources. 

The  receiver  depth  is  not  included  in  the  system  of 
equations,  and  any  number  of  receivers  can  therefore  he  treated 
with  only  one  solution. 
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MODELLING  OP  PULSE  PROPAGATION  IN  LAYERED  MEDIA  USING 
A  NEW  FAST  FIELD  PROGRAM 
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19026  La  Spezla,  Italy 


ABSTRACT 

Numerical  modelling  of  acoustic  pulse  propagation  In  stra¬ 
tified  solid  media  Is  often  found  to  be  impractical  due  to  the 
comprehensive  calculations  involved.  Here  a  new  numerical  model, 
of  the  fast  field  type.  Is  presented.  Instead,  of  using  the 
Thomson-Haskell  matrix  method,  as  done  In  earlier  models  of  the 
same  type,  the  depth-separated  wave  equation  Is  solved  by  a 
numerical  technique  very  similar  to . that  used  In  finite  element 
programs.  The  speed  Improvement  has,  been  considerable,  especially 
in  cases  with  many  sources  and  receivers.  The  model  has  been 
implemented  an  an  FPS164  array  processor  and  used  for  analysis  of 
seismic  pulse  propagation  in  a  shallow  water  environment  and 
reflection  of  pulsed  ultrasonic  beams  from  a  fluid-solid  inter¬ 
face. 


INTRODUCTION 

A  number  of  numerical  models  are  available  for  investigating 
sound  propagation  in  layered  media.  Each  model  is  based  on  a  set 
of  assumptions  and  approximations  and  dedicated  to  specific  appli¬ 
cations.  In  relation  to  underwater  acoustics  four  types  of  models 
are  of  interest,  known  as  ray,  normal-mode,  parabolic  equation, 
and  fast  field  models.  Jensen  <1>  has  reviewed  and  classified 
these  models  and  here  we  will  concentrate  on  models  of  the  last 
type,  the  fast  field  programs  (FFP) . 

These  models  yield  an  exact  solution  to  the  depth-separated 
wave  equation  for  horizontally  stratified  environments.  Depending 
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6  SUMMARY  AND  CONCLUSIONS 

We  have  briefly  presented  an  overview  over  the  moat  commonly 
used  propagation  models  In  underwater  acoustics  pointed  out  their 
areas  of  applicability,  and  demonstrated  their  ability  to  accura¬ 
tely  describe  acoustic  propagation  in  complicated  ocean  environ¬ 
ments.  It  has  also  been  shown  that  a  good  agreement  between 
theory  and  experimental  data  can  be  obtained  only  by  including 
such  features  as  bottom  layering,  bottom  rigidity  (shear),  scat¬ 
tering  at  rough  boundaries,  range-varying  environments,  etc. 
Hence  the  full  complexity  of  a  real  ocean  environment  must  be  con- 
sidered  In  the  numerical  models  for  accurately  predicting  the  pro¬ 
pagation  conditions  In  a  given  area  for  a  broad  range  of  source 
frequencies.  Finally,  the  general  applicability  of  the  numerical 
models  has  been  demonstrated  by  applying  the  models  to  some  basic 
wave-propagation  problems  not  tractable  by  analytical  methods. 
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Fig.  20  Down-slope  propagation  for: 

a)  a  5°  slope  followed  by,  a  flat  bottom,  and 

b)  a  smoothly  changing  bottom  slope. 

There  are  several  ways  to  Increase  the  degree  of  coupling 
between  modes  for  down-slope  situations.  Besides  the  Increased 
coupling  with  frequency  and  slope  angle,  coupling  can  also  be 
caused  by  abrupt  changes  In  bottom  slope.  This  is  shown  in 
Fig.  20a,  which  Is  the  same  environment  as  In  Fig.  19a,  except 
that  the  bottom  now  Is  flat  beyond  5  km.  The  abrupt  change  In 
bottom  slope  at  5  km  results  in  a  complicated  contour  pattern 
indicating  interference  between  modes  and,  hence,  the  initial 
first  mode  has  generated  (or  coupled  into)  higher  modes  after  the 
vertex.  In  this  case  as  many  as  six  modes  can  exist  in  the  flat 
part  beyond  5  km.  That  mode  coupling  is  associated  with  the  sharp 
change  of  the  vertex  is  shown  in  Fig.  20b  where  we  gradually  go 
from  a  5°  slope  to  a  flat  bottom.  We  see  that  the  propagating 
mode  now  adapts  to  the  changing  water  depth  exhibiting  an 
"adiabatic''  behavior. 
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Fig.  19  Down-slope  propagation  over  a  constant  5°  slope 

a)  25  Hz,  no  mode  coupling, 

b)  500  Hz,  strong  mode  coupling . 

Fig.  19a  shows  that  some  energy  propagates  straight  into  the 
bottom  at  short  ranges  (coupling  into  the  continuous  spectrum). 
However,  beyond  the  nearfield,  propagation  within  the  wedge  is 
clearly  adiabatic  with  the  one  propagating  mode  adapting  well  to 
the  changing  water  depth.  At  range  20.  km  the  energy  is  entirely 
contained  in  the  local  first  mode,  even  though  as  many  as  21  modes 
can  exist  in  a  water  depth  of  1800  m. 

By  increasing  the  frequency  to  500  Hz  (Fig.  19b),  strong  mode 
coupling  occurs  within  the  wedge.  At  long  ranges,  two  inter¬ 
ference  nulls  are  present  in  the  energy  distribution  over  depth 
indicating  that  the  energy  is  now  partitioned  among  a  few  lower- 
order  modes. 
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method*  However,  his  approach  allows  for  only  one  source/recei ver 
combination  at  a  time. 

In  the  SAFARI  model  Eq.  I?  is  set  up  and  solved  directly. 
When  the  layer  series  is  given,  it  is  possible  to  determine  a 
priori  a  mapping  between  the  equations  to  be  satisfied  at  each 
interface  and  the  global  system  of  equations,  Eq.  19.  The  unknown 
arbitrary  functions  are  mapped  into  the  vector  A^ ,  the  coefficients 
into  the  matrix  Cj j ,  and  the  source  contributions  into  the  right 
side  Rj  oi  Eq.  19.  If  more  than  one  source  is  present,  the  contri¬ 
butions  are  simply  added. 

The  receiver  depth  is  not  included  in  Eq.  19,  thus  yielding 
the  possibility  of  determining  the  field  at  any  depth  from  Eqs.  12 
to  16  with  only  one  solution  of  Eq.  19. 

The  mapping  between  the  local  and  the  global  system  of 
equations  is  very  similar  to  the  technique  used  in  finite  element 
programs.  By  using  this  technique  the  computer  code  can  become 
very  efficient  since  a  significant  number  of  calculations  can  be 
done  only  once.  Furthermore,  the  code  will  be  straightforward  to 
vectorize,  thus  making  it  well  suited  for  implementation  on  an 
array  processor. 

The  solution  of  Eq.  19  is  performed  by  means  of  gaussian  eli¬ 
mination  with  partial  pivoting. 

In  cases  with  only  one  source/receiver  combination,  SAFARI  is 
found  to  be  an  order  of  magnitude  faster  'than  Kutschale 's  model. 
In  cases  with  several  sources  or  receivers  the  speed  gain  is  of 
course  much  bigger. 

This  surprisingly  high  speed  gain  is  believed  to  be  due  to  the 
mapping  technique,  which  partly  ensures  that  a  large  amount  of 
unnecessary  calculation  is  avoided,  and  partly  makes  efficient 
programming  possible. 

The  implementation  of  SAFARI  on  an  FPS164  attached  processor 
has  given  at  least  one  order  of  magnitude  further,  and  a  number  of 
problems,  which  were  impractical  to  treat  numerically,  can  now  be 
solved  with  an  acceptable  amount  of  computation  time.  A  couple  of 
examples  are  given  in  the  following. 


3  MODELLING  OF  SEISMIC  INTERFACE  WAVES 

The  Importance  of  the  shear  properties  of  the  sea-bed  for 
acoustic  wave  propagation  in  shallow  water,  is  well  established  <5>. 
Unfortunately  the  shear  parameters  are  very  difficult  to  Isolate 
experimentally.  The  shear  parameters  are,  however,  indirectly  pre- 
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Fig.  1  Stratified  shallow  water  environment 


sent  through  the  properties  of  the  measurable  seismic  interface 
waves,  and  much  effort  has  been  put  into  an  experimental  investiga¬ 
tion  of  these  <9>,  <10>. 

Since  no  applicable  inverse  models  are  available  the  shear 
properties  are  determined  by  "trial-and-error”  —  methods  using 
numerical  propagation  models  <10>,  <11>.  Usually  several  parame¬ 
ters  are  unknown  and,  since  the  calculations  needed  for  each  com¬ 
bination  are  rather  comprehensive,  $he  determination  of  the  shear 
parameters  in  this  way  can  become  very  expensive  in  terms  of  calcu¬ 
lation  time.  It  is  therefore  important  a  priori  to  be  able  to 
determine  approximate  values  directly  from  the  experimental  data. 

With  this  in  view,  the  SAFARI  model  has  been  used  to  clarify 
some  of  the  propagation  characteristics  of  the  seismic  interface 
waves  for  different  shallow  water  epvironments.  A  detailed 
description  of  the  investigation  is  given  In  <I2>,  and  only  a 
couple  of  examples  will  be  given  here. 

In  order  not  to  obscure  the  basic  principles,  a  simple 
2-layered  model  was  chosen  for  the  sea-bed.  Fig.  1.  Below  the 
water  column  of  depth  c^,  a  single  sediment  layer  of  thickness  d8 
covers  a  half-space  of  rock  or  rock-like  material.  The  bottom 
materials  used  in  the  examples  and  their  assumed  properties  are 
listed  in  Table  1.  The  complex  Lame  constants  are  not  given  expli¬ 
citly  in  the  table.  Instead  the  compressions!  and  shear  velocities 
are  shown,  together  with  their  respective  dampings  in  dB  per  wave¬ 
length. 
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TABLE  1 


MATERIAL  PROPERTIES 


Mat  •rial 

Density 

P(|/CB*J 

Comp re s  s i onaj 
Speed 
ct(m/fc) 

Shea  r 
Speed 
CT(m/s) 

Compress  lona  1 
attenuat i on 

Shea  r 

attenuat ion 

Vt<«ib/a> 

Water 

1.0 

1500 

- 

- 

Silt 

IS 

1600 

200 

1.0 

2  0 

Sand 

2  0 

1800 

600 

0.7 

l  5 

Limestone 

2  2 

2250 

1000 

0.4 

1  Q 

Basalt 

2.6 

5250 

2500 

0.2 

0  5 

The  water  depth  Is  chosen  to  be  100  m,  and  a  silt  layer  of 
50  m  thickness  is  combined  with  either  a  limestone  or  a  basalt  sub¬ 
bottom. 

A  point  source  is  placed  in  the  middle  of  the  water  column  at 
50  m  depth,  and  the  vertical  particle  velocity  at  the  top  of  the 
sediment  layer  is  calculated  using  a  pressure  amplitude  of  1  Pa  at 
a  distance  of  1  m  from  the  source.  Only  frequencies  below  and 
around  the  cut-off  frequency  for  the  first  water  mode  are  con¬ 
sidered,  which  in  the  present  case  means  frequencies  below  10  Hz. 
In  the  first  test  case  a  50  m  silt  sediment  layer  covers  a 
limestone  sub-bottom.  Figure  2  shows  the  modulus  of  the  integrand 
in  the  Hankel  transform  of  the  vertical  velocity  at  the  top  of  the 
sediment  layer  for  a  frequency  of  1.5  Hz.  No  water  modes  are  pre¬ 
sent  at  this  frequency,  but  two  mode-like  peaks  can  be  observed  at 
wavenumbers  corresponding  to  phase  velocities  of  300  m/s  (peak  '1') 
and  840m/s  (peak  '2').  These  peaks  correspond  to  interface  waves, 
and  will  be  denoted  the  first  and  second  interface  mode,  respec¬ 
tively.  The  first  Interface  mode  is  best  excited  (highest 
amplitude),  but  it  also  has  the  highest  damping  (widest  peak). 
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Fig.  3  Vertical  particle  velocity  as  function  of  range  at  1.5  Hz. 
50  m  silt  on  limestone. 


Figure  3  shows  the  corresponding  vertical  particle  velocities  at 
ranges  up  to  50  km.  Due  to  the  high  damping  of  the  first  interface 
mode,  its  contribution  is  significant  only  for  ranges  shorter  than 
2  km,  beyond  which  the  second  Interface  mode  becomes  dominant. 


If  the  source  is  not  of  stationary  type,  but  transient,  the 
different  velocities  of  the  two  interface  modes  will  yield  dif¬ 
ferent  arrival  times,  and  the  presence  of  the  first  interface  mode 
will  be  measurable  also  at  larger  ranges.  The  phase  and  group 
velocities  of  the  two  Interface  modes  have  been  calculated  in  the 
frequency  band  of  Interest,  0.1  to  10  Hz,  with  a  resolution  of 
0.1  Hz.  The  results  are  shown  in  Fig.  4  together  with  an  excita¬ 
tion  measure,  which  somewhat  arbitrarily  has  been  chosen  to  be  the 
particle  velocity  at  10  km  range.  » 
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Fig.  4  Excitation  and  dispersion  curves.  50  m  silt  on  limestone. 
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There  is  only  one  Interface  mode  at  frequencies  below  1  Hz. 
It  is  slightly  dispersive,  with  phase  and  group  velocities 
approaching  those  of  a  Rayleigh  wave  on  a  limestone  half-space.  At 
1  Hz  a  very  sharp  transition  zone  appears,  where  the  velocities 
drop  dramatically  to  values  approaching  those  of  a  Scholte  wave  at 
a  water /silt  interface.  At  1.8  Hz  the  group  velocity  reaches  a 
minimum  of  100  m/s,  i.e.  half  the  shear  speed  in  silt.  The  second 
interface  mode  has  a  sharp  cut-off  at  the  transition  frequency,  and 
afttrr  a  distinct  minimum  of  the  group  velocity  it  appears  as  a 
logical  continuation  of  the  low-frequency  part  of  the  first  inter¬ 
face  mode.  Above  2  Hz  the  excitation  (solid  line)  of  the  second 
interface  mode  decreases  due  to  the  increased  distance  in  terms  of 
wavelengths  of  the  source  from  the  silt/limestone  interface,  along 
which  the  second  interface  mode  propagates.  The  sharp  peak  on  the 
excitation  curve  around  4.5  Hz  is  the  first  propagating  water 
mode. 


The  existence  of  mode  transition  zones  is  well  known  from  the 
theory  of  vibration  of  elastic  plates,  Mindlin  <13>,  where  they 
appear  near  the  thickness-shear  frequencies.  These  are  the  fre¬ 
quencies  at  which  an  infinite  elastic  plate  can  perform  free  shear 
vibrations  with  vanishing  vertical  displacements.  Now  consider  the 
silt  layer  as  an  infinite  elastic  plate  in  welded  contact  with  an 
infinite  rigid  half-space.  The  first  thickness-shear  frequency 
would  then  correspond  to  that  of  a  free  silt  plate  of  the  double 
thickness,  Mindlin  <13>: 


(22) 


where  Gy  is  the  shear  velocity  and  dg  Is  the  thickness  of  the 
layer.  If  the  parameters  for  silt  are  applied  to  Eq.  22,  we  obtain 
f-pg  equal  to  1  Hz,  which  is  very  close  to  the  observed  transition 
frequency  in  Fig.  4. 


Since  the  transition  frequency  can  be  observed  in  experimental 
results,  its  correlation  with  the  thickness-shear  frequency  could 
yield  a  direct  method  of  determining  the  shear-wave  velocity  in  a 
single  sediment  layer  overlying  a  rigid  half-space. 


To  summarize  the  general  propagation  characteristics  for  this 
test  case  (Fig.  4),  the  computed  particle  velocity  at  10  km  range 
is  entirely  associated  with  the  first  interface  mode  below  the 
transition  frequency  (1  Hz).  This  Interface  mode  is  strongly 
related  to  a  Rayleigh  wave  on  a  limestone  half-space  below  1  Hz, 
while  it  becomes  an  interface  wave  connected  with  the  water/silt 
interface  at  frequencies  above  1  Hz.  Zn  this  frequency  regime  the 
second  interface  mode  appears,  and  it  is  mainly  related  to  the 
sllt/llmestone  interface.  Finally,  the  first  water-borne  mode 
appears  at  around  4.5  Hz. 
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Fig.  5  Excitation  and  dispersion  curves.  50  m  silt  on  basalt. 


A  silt  layer  an  basalt  was  also  studied  in  order  to  analyze 
the  effect  of  different  sub-bottom  materials.  In  this  case 
(Fig.  5)  the  low-frequency  part  of  the  first  interface  mode  is 
excited  only  slightly,  but  the  transition  zone  can  again  be 
observed  near  the  thickness- shear  frequency  of  1  Hz.  Above  1  Hz 
the  dispersion  curves  are  very  similar  to  those  obtained  earlier 
for  a  sub-bottom  of  limestone  (Fig.  4). 

Synthetic  seismograms  have  been  produced  for  the  two  test 
cases  in  order  to  illustrate  the  time-rdomain  effect  of  the  features 
described  above.  The  source  is  assumed  to  be  half  a  sine  wave  of 
1.5  Hz  sent  through  an  ideal  0.5  to *3.2  Hz  band-pass  filter.  This 
frequency  range  has  been  chosen  as  it  contains  frequencies  on  both 
sides  of  the  transition  frequency  for  a  silt-sediment  layer  of  50  m 
thickness.  The  transfer  functions  were  calculated  to  a  resolution 
of  0.01  Hz  and  multiplied  by  the  spectrum  of  the  source.  The  time 
series  were  then  created  by  means  of  the  fast-Fourier  transform  at 
ranges  of  1,  2,  3,  4  and  5  km. 

Figure  6  shows  the  result  for  test  case  1,  i.e.  a  50  m  silt 
layer  on  a  limestone  half-space.  The  first  weak  arrival 
corresponds  to  a  compressional  wave  in  the  limestone,  whereas  the 
first  significant  arrival  corresponds  to  the  Rayleigh  wave  velocity 
of  the  limestone  (900  m/s).  This  arrival  consists  of  the  low- 
frequency  parts  of  the  first  and  second  interface  modes.  A  clear 
dispersion  can  be  observed  corresponding  to  the  distinct  minimum  in 
the  group  velocity  of  the  second  interface  mode  in  Fig.  4.  The 
slow  highly  damped  wave  corresponding  to  the  first  interface  mode 
above  the  transition  frequency  propagates  with  group  velocities 
between  100  and  180  ra/s,  again  in  good  agreement  with  Fig.  4.  At 
ranges  greater  than  4  km  this  arrival  has  negligible  amplitude. 
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Fig.  6 

Stacked  synthetic  seismograms. 
50  m  silt  on  limestone. 
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Fig.  7 

Stacked  synthetic  seismograms . 
50  m  'Silt  on  basalt. 


As  described  above,  the  low-frequency  part  of  the  first  Inter¬ 
face  mode  will  not  be  significantly  excited  If  the  sub-bottom  Is 
basalt.  The  maximum  excitation  at  10  km  range  lies  at  2.5  Hz 
(Fig.  5)  and  Is  due  to  a  strong  excitation  of  the  second  mode. 
These  properties  are  also  reflected  in  the  synthetic  seismograms 
(Fig.  7).  The  fastest  arrival,  apart  from  the  weak  compresslonal 
wave,  has  its  major  frequency  content  above  2  Hz  and  corresponds 
mainly  to  the  high-frequency  end  of  the  second  Interface  mode.  The 
severe  dispersion  of  the  corresponding  arrival  in  Fig.  6  Is  not 
present  here,  and  the  slow  part  of  the  first  Interface  mode  is  well 
separated,  travelling  at  group  velocities  between  100  and  200  m/s. 
The  synthetic  seismograms  in  Figs.  6  and  7  are,  at  least  qualitati¬ 
vely,  very  similar  to  those  observed  during  experiments.  Fig.  8 
<15>. 

With  4096  sample  points  in  the  wavenumber  space  the  calcula¬ 
tion  time  on  an  FPS164  array  processor,  was  9  secnds  for  each  fre¬ 
quency  or  40  minutes  in  total  for  the  synthetic  seismograms  in 

Figs.  6  or  7. 
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Fig.  8  Stacked  experimental  seismograms ,  <1 5>. 


4  REFLECTION  OF  A  PULSED  BEAM  AT  A  FLUTD/ SOLID  INTERFACE 

As  mentioned  above,  the  solution  technique  used  in  SAFARI 
yields  the  possibility  of  treating  problems  in  which  many  sources 
and  receivers  are  involved.  , 

Jensen  <1>,  used  the  model  to  treat  the  problem  of  reflection 
of  narrow  ultrasonic  beams  at  a  water/aluminium-oxide  interface  at 
grazing  angles  near  the  Rayleigh  angle,  i.e.  the  angle  where  a 
leaky  interface  wave  Is  excited.  A  monochromatic  analysis  was  per¬ 
formed  in  order  to  clarify  the  influence  of  beam  divergence  on  the 
reflection  pattern,  and  excellent  agreement  was  obtained  with  both 
experimental  and  theoretical  results  reported  in  the  literature. 
In  addition,  some  apparent  discrepancies  between  experimental  and 
theoretical  results  were  resolved. 

One  of  the  results  from  <l>  is  shown  in  Fig.  9.  A  parallel 
gaussian  beam  is  generated  by  a  linear  vertical  source-array  of  649 
elements,  spaced  half  a  wavelength  apart.  The  mid-point  of  the 
array  is  placed  40  cm  above  the  interface,  the  frequency  is  2  MHz, 
and  the  grazing  angle  of  the  beam  is  75.2”,  corresponding  to  the 
Rayleigh  angle  for  water /aluminium-oxide.  The  splitting  of  the 
reflected  beam  is  easily  observed.  The  specularly  reflected  beam 
and  the  beam  caused  by  the  leaky  interface  wave  are  separated  by  a 
strip  of  low  intensity. 
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Fig.  9 

Reflection  of  a 
2  MHz  beam  from 
a  water /aluminum- 
oxide  interface 
at  the  Rayleigh 
angle. 


As  pointed  out  by  Bertoni  and  Tamir  <14>  this  zero-strip  is 
due  to  destructive  interference  between  the  two  parts  of  the 
reflected  field.  If  the  surfaces  of  equal  phase  are  plane  and  per¬ 
pendicular  to  the  direction  of  propagation,  which  is  the  case  for 
parallel  beams,  the  two  reflected  beams  thus  have  to  be  180°  out  of 
phase. 


This  feature  can  of  course  be  demonstrated  by  making  a  cut 
perpendicular  to  the  propagation  direction,  but  if  the  continuous 
beam  is  replaced  by  a  pulsed  beam,  the  phase  shift  will  yield  an 
amplitude  inversion  of  the  pulse  in  the  two  beams,  Independent  of 
the  position  of  the  receivers  within  the  beams.  Furthermore,  a 
pulse  calculation  will  yield  information  on  arrival  times  for  the 
pulse  at  different  receivers. 

To  demonstrate  this  the  pulse  version  of  SAFARI  was  used  to 
calculate  the  received  signals  at  five  different  receivers,  all 
situated  20  cm  above  the  interface  as  indicated  in  Fig.  9. 

As  pulse  calculations  in  SAFARI  are  performed  by  means  of 
discrete  Fourier  synthesis,  it  is  very  important  to  choose  the  time 
window  correctly  in  order  to  avoid  time-domain  aliasing. 

The  approximate  arrival  times  were  therefore  determined  using 
a  narrow-banded  source  pulse  (1.8  to  2.2  MHz).  The  time  window  was 
chosen  to  be  0  to  1  ms,  yielding  frequency  steps  of  1  kHz.  The 
calculated  time  series  for  the  five  receivers  are  shown  in  Fig.  10. 
The  vertical  scale  is  arbitrary,  but  is  the  same  for  all  five 
receivers. 
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Fig.  lO 

Received  narrow-band  pulses 
at  receivers  I-V. 


The  pulse  arrives  at  the  reference  receiver  (I)  at  t  ~  140  ps. 
As  could  be  expected  from  Fig.  9,  no  significant  arrivals  appear  at 
receiver  II,  whereas  the  reflected  pulses  arrive  at  the  three  last 
receivers  at  times  between  4 10  and  440  ps  approximately. 

A  closer  analysis  Is  now  possible  using  time  windows  of  125  us 
width.  In  order  to  obtain  short  pulses,  the  bandwidth  is  Increased 
to  1.5  MHz  to  3.0  MHz.  The  results  are  shown  for  receiver  I,  III 
and  IV  In  Fig.  11.  The  time  scale  ’shows  the  time  relative  to  the 
expected  arrival  times.  These  are  determined  by  geometrical  means. 
A  plane  wave  of  grazing  angle  75.2°  that  passes  the  reference 
receiver  I  at  t  ■  139.30  ps  will  be  reflected  and  pass  receiver  III 
at  t  *  416.00  ps.  The  calculated  pulse  at  this  receiver  Is  seen  to 
correspond  exactly  to  this  behaviour;  the  specularly  reflected  beam 
has  no  phase  shift.  The  loss  In  amplitude  is  due  to  the  fact  that 
same  of  the  energy  is  transferred  int6  the  leaky  Rayleigh  wave. 
The  expected  arrival  at  receiver  IV  is  determined  In  a  slightly 
different  way.  A  ray  at  a  grazing  angle  of  75.2°  Is  assumed  to 
pass  the  reference  receiver  at  t  ■  139.30  ps.  When  It  reaches  the 
interface  it  travels  at  the  speed  of  the  Rayleigh  wave 
(5825.55  m/a)  along  the  interface,  and  finally  it  travels  upwards 
at  angle  of  75. 2#  to  receiver  IV.  The  arrival  time  determined  by 
these  assumptions  is  424.60  ps,  again  in  excellent  agreement  with 
the  numerical  result  in  Fig.  11.  The  180°  phase  shift  is  easily 
observed,  thus  verifying  the  explanation  for  the  2ero-strip  given 
in  <14>. 
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Fig.  11  Arrival  of  broad-band  pulses  at  receivers  I,  III  and  IV 
related  to  expected  arrival  time  tQ. 

With  512  sampling  points  in  the  wavenumber  space,  the  calcula¬ 
tion  time  on  the  FPS164  was  8  seconds  for  each  frequency,  in  total 
40  minutes  for  Figs.  10  and  25  minutes  for  Fig.  11. 


CONCLUSIONS 

* 

A  new  fast  field  program,  SAFARI,  has  been  developed.  The 
program  is  meant  for  modelling  sound  and  stress  wave  propagation  in 
horizontally  stratified  fluid  and  solid  media.  By  introducing  a 
more  efficient  solution  technique,  the  computational  speed  has  been 
Improved  by  an  order  of  magnitude,  in  some  cases  even  more,  com¬ 
pared  with  earlier  models  of  the  same  type.  Furthermore,  the 
SAFARI  model  is  capable  of  treating  several  'sources  and  receivers 
with  only  one  solution  of  the  wave  equation,  thus  making  it 
feasible  to  treat  problems,  in  which  the  field  is  generated  by  ver¬ 
tical  source  arrays. 

The  model  is  basically  monochromatic,  but  the  increase  in 
speed,  together  with  the  fact  that  the  model  is  well  suited  for 
implementation  an  an  array  processor,  has  yielded  the  possibility 
of  treating  pulse  propagation  by  means  of  Fourier  synthesis  with 
limited  demands  on  computation  time.  This  has  been  demonstrated  by 
a  couple  of  examples,  but  the  range  of  problems  that  could  be 
treated  is  of  course  much  wider. 
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